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SYMBOLS 

Sections  II  and  III 
a  local  speed  of  sound 

a8  critical  speed  of  sound 

M  Maeh  Number 

M8  critical  Mach  number 

q  velocity  magnitude 

rc  initial  radius  of  curvature  ©f  expanded  free  jet  boundary 

R  ratio  of  distance  from  center  of  expansion  t©  distance  ©f 

center  of  axis  ©f  symmetry 
u,v  velocity  components 

y  Cartesian  coordinates 

constants  ©f  integration 

Ki,$2 

Y  specific  heat  rati© 

<j>  polar  coordinate 

A  C(Y-l)/(Y+!))% 

p  local  Mach  angle 

p  auxiliary  angle 

A  angle  at  the  center  containing  the  region  of  centered 

expansion 

5  angle  ©f  deflection  of  the  streamline  at  the  center  of 

expansion 

9  streamline  angle 

^  4>-B ! 


2 


SYMBOLS 


M 

M» 

R 

Rw 

# 

u?v 
X,  Y 


z,r 


n 


c 


F 

M 


Section  IV 

sonic  velocity 
Mach  number 
critical  Mach  number 

normalized  throat  radius  of  curvature,  R=Rw/rs 

nozzle  throat  radius  of  curvature 

nozzle  throat  radius 

velocity  components 

Cartesian  coordinates 

normalized  coordinate  of  Hall's  expansion 

Dutton-Addy  expansion  variable 
streamline  angle 


Subscripts 

point  in  nozzle  where  the  circular  throat  joins  the 
conical  section 

conditions  at  final  expansion  fan  line  as  R-»-0 
model 
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F 

L 

M 

M 

P 


Subscripts 

conditions  at  final  expansion  fan  line  as  R-+0 
conditions  at  initial  expansion  fan  line  as  R-H) 
model 

limiting  value  as  the  center  of  the  expansion  is  approached 

along  the  nozzle  boundary  streamline 

prototype 
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I.  INTRODUCTION 


The  interaction  of  a  rocket's  plume  with  the  external  flow  can  signifi¬ 
cantly  affect  the  aerodynamic  performance  of  tre  rocket.  The  adverse  effects 
are  well  known  and  include  increased  drag,  base  heating  and  plume  induced 
separation.  There  is  at  least  one  instance  whi  re  plume  interaction  can  be 
beneficial.  Some  rocket  configurations  induce  large  stable  pitching  moments 
in  the  transonic  flight  regime  which  cause  the  vehicle  to  be  overly  wind  sen¬ 
sitive.  The  interaction  between  the  slipstream  and  an  underexpanded  plume  can 
act  to  de- stabilize  the  rocket  in  this  region. 

Whether  to  avoid  the  adverse  conditions  or  to  take  advantage  of  beneficial 
effects,  the  missile  designer  needs  methods  of  investigating  the  plume- 
slipstream  interaction.  In  the  transonic  region,  experimental  testing  is  the 
primary  method  of  obtaining  data  and  it  is  most  desirable  to  have  wind  tunnel 
models  which  simulate  the  proper  plume  characteristics  without  actual  pro¬ 
pellant  burning.  Korst^-  has  presented  a  method  for  designing  models  of  coni¬ 
cal  nozzles  which  produce  the  same  geometrical  plume  as  the  prototype  nozzles 
while  accounting  for  significant  viscid  and  inviscid  aspects  of  the  base  flow 
problem.  The  flow  analysis  used  by  Korst  is  based  on  concepts  developed  by 
Johannesen  and  Meyer^  which  allow  the  flow  field  near  the  centered  expansion 
at  the  nozzle's  exit  to  be  expressed  in  the  form  of  a  series  expansion  with 
respect  to  the  radius  vector.  For  a  set  of  specified  prototype  nozzle  con¬ 
ditions,  Korst 's  modeling  theory  can  be  used  to  calculate  model  nozzles  with 
the  proper  combination  of  exit  lip  Mach  number,  wall  flow  acceleration  and  lip 
wall  angle.  For  nozzles  with  nominal  divergent  sections  where  source  flow  is 
approximated  the  modeling  requirements  are  reduced  to  determining  the  exit 
Mach  number  and  exit  wall  angle.  This  modeling  concept  has  also  been  pre¬ 
sented  in  a  paper  by  Korst  and  Deep^  which  suggests  a  numerical  solution  for 
one  of  the  approaching  flow  velocity  terms  as  an  alternative  to  the  closed 
form  solution  of  Johannesen  and  Meyer. 

This  latter  scheme,  while  determining  the  requisite  exit  conditions,  does 
not  account  for  the  details  of  the  nozzle  geometry  necessary  to  produce  these 
exit  conditions.  For  a  conical  nozzle  with  a  circular  arc  throat  section  the 
geometrical  quantities  which  must  be  determined  are  throat  radius  of  cur¬ 
vature,  nozzle  length  and  the  axial  position  where  the  circular  arc  section 
and  conical  divergent  section  join.  To  evaluate  the  effects  of  throat  radius 
of  curvature  use  is  made  of  the  transonic  throat  flow  approximations  origi¬ 
nally  proposed  by  Hall^  and  corrected  by  Kliegel  and  Levine^  and  later  refined 
by  Dutton  and  Addy^.  The  flow  field  as  calculated  by  this  method  then  provi¬ 
des  starting  conditions  for  a  method  of  characteristics  routine  which  solves 
the  flow  field  in  a  radial  plane  fashion  until  the  exit  conditions  specified 
by  Korst 's  modeling  theory  are  matched.  The  point  where  the  exit  conditions 
match  is  then  the  nozzle  length.  The  point  where  the  circular  arc  section  and 
conical  section  join  can  be  determined  from  purely  geometrical  considerations. 

Three  programs  written  by  Prof.  Korst  in  BASIC  for  the  Hewlett-Packard 
9830  were  re-written  in  FORTRAN  by  the  author  and  combined  into  a  single 
interactive  model  design  program.  This  report  presents  an  overview  of  the 
Johannsen  and  Meyer  theory,  the  throat  flow  approximation  theory  and  Korst' s 
modeling  theory.  The  report  also  illustrates  how  the  three  theories  are  used 
to  develop  the  model  design  program.  The  program  also  calculates  plume  shape 
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using  a  method  of  characteristics  routine  which  is  initialized  by  the  exit 
conditions  determined  from  the  modeling,,  The  FORTRAN  program  is  included  with 
user  instructions  and  a  sample  case. 

II,  THEORY  OF  JORANNESEN  AND  MEYER 

The  Korst  plume  modeling  theory  is  based  on  the  theory  of  Johannesen  and 
Meyer-  which  provides  an  approximate  method  for  solving  the  flow  field  at  the 
lip  of  an  axially-symmetric  nozzle.  The  approach  assumes  the  flow  is  a  per¬ 
fect  gas  that  is  isentropic,  irrotational  and  steady. 

When  a  uniform,  two-dimensional  flow  with  these  properties  expands  around 
a  corner  it  is  called  Prandtl-Meyer  flow  and  is  characterized  by  generating 
centered  simple  waves.  In  simple-wave  flow  the  disturbances  created  by  the 
corner  are  propagated  as  pressure  waves  along  the  Mach  lines.  Further,  the 
Mach  lines  are  assumed  straight  and  the  flow  properties  along  the  Mach  lines 
are  uniform.  The  velocity  magnitude  at  any  point  in  a  centered  wave  is 
characterized  by  streamlines  such  that  all  the  Mach  lines  pass  through  a  com¬ 
mon  origin.  When  the  flow  is  axially-symmetric,  the  expansion  is  no  longer 
simple,  though  it  is  still  centered. 

The  method  proposed  by  Johannesen  and  Meyer  solves  the  flow  field  in  the 
neighborhood  of  the  lip  by  assuming  that  the  flow  may  be  divided  into  three 
regions,  a  centered  expansion  occurring  at  the  lip,  and  flow  fields  before  and 
after  the  expansion,  A  polar  coordinate  system  with  variables  R  and  <j)  is 
employed,  see  Figure  1,  where  R  is  the  ratio  of  the  distance  from  the  center 
of  the  expansion  to  the  distance  of  the  center  of  the  axis  of  symmetry,  such 
that  the  exit  radius  of  the  nozzle  has  a  value  of  1.  The  polar  coordinate,  <S>3 
is  measured  from  the  exit  plane  of  the  nozzle,  A  series  expansion  in  powers 
of  R  gives  the  following  velocity  components : 

u(R, <j>)  =  u0  O)  +  Rup  (4>)  +  R2U2  (<t>)  +  ...  (1) 

v(R,  <j>)  =  v0  (<J>)  +  Rvp  (4>)  +  R2V2  (<t>)  +  ...  (2) 

The  first  term  in  this  expansion  gives  a  solution  identical  to  the  two- 
dimensional  solution,  that  is  a  uniform  flow  field  before  and  after  the  expan¬ 
sion  and  a  Prandtl-Meyer  flow  field  in  the  region  of  the  expansion.  The 
equations  for  the  second  term  in  the  series  must  be  integrated  to  provide  a 
solution.  With  the  resulting  equations  and  the  proper  application  of  boundary 
conditions,  the  initial  curvature  of  a  jet  boundary  can  be  calculated  with  an 
error  of  0(R2). 

With  the  further  assumption  that  the  expansion  at  the  lip  is  shock-free, 
the  development  of  Johannesen  and  Meyer  is  presented  in  the  following 
paragraphs. 

The  polar  equations  of  continuity  and  of  motion  may  be  combined  to  give 
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where  a  is  local  speed  of  sound,  a*  is  critical  speed  of  sound,  y  is  specific 
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heat  ratio  and.  \  »^Jy+1 
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The  irrotationality  of  the  flow  is  expressed  by; 


R9cf> 


1  3(Rv) 

R  9R 


(5) 


If  the  series  representation  of  the  velocity  components  as  given  by  (1)  and 
(2)  are  substituted  into  equations  (4)  and  (5),  coefficients  of  various  powers 
of  R  are  equated  and  O(R^)  and  greater  are  discarded  the  following  equations 
are  obtained; 
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where  u'  =  — rr  and  ul*  =  — rr 
o  acp  i  dtp 

If  equation  (3)  Is  expanded  with  the  series  representations  of  u  and  v,  0(r2) 
and  greater  are  again  discarded,  powers  of  R  are  collected  and  equations  (6) 
to  (9)  are  used  the  following  two  equations  can  be  derived; 
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From  equation  (10)  it  is  seen  that  there  are  two  solutions  to  the  system  of 
equations.  One  solution  approaches  uniform  flow  and  the  other  solution 
approaches  a  Prandtl-Meyer  expansion  as  R-»0, 


A,  The  Centered  Expansion, 
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(12) 


To  arrive  at  equations  for  the  centered  expansion  let 


2  2 

v  -a  =0 
o  o 

then  from  (6)  and  (8)  after  integrating 

u0(<t>)  =  X~  sin  [x  (4.  +  pQ)]  (13) 

v0(4>)  =  a*  cos  [x  ( 4>  +  f30)]  (14) 

where  the  constant,  P0,  will  be  determined  by  satisfying  the  initial  con¬ 
ditions  of  the  approaching  flow. 

Equation  (11)  can  be  reduced  by  substituting  in  equation  (12)  after  re¬ 
arranging  equation  (11)  becomes: 

2 


2u  u.  + 
o  1 


["  3y-l  /uo^ 

-  1 

Y+l  \V  / 

L.  i  O 

v  u.  -  v  (v  sin<j>  -  u  cost}))  -  0 
o  1  o  o  o 


(15) 


Then  after  substituting  (13)  and  (14)  into  (15)  and  integrating,  Johannesen  and 
Meyer  derived  the  following: 


u.  (<j>)  =  -a*  (cosn) 
1  2  X 
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where  q  =  X  (4>  +  p0)  and  is  an  arbitrary  constant  of  integration,  and 

J1  -V  X  2 

I-^d)  =J  cos  (y/X)  (siny)  2  (cosy)  2  dy 

ni 

2 

I2  (n)  =  f  n  sin  (y/X)  (siny)-1  (cosy)  ^  Y  dy 

n.- 


de) 


-3/2 


f-3 

be 
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I  a  (in)  =f11  cos  (y/X)  (sin  y)'3/(2cosy)"2Y  "  d7 
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where  y  is  the  Cartesian  coordinate  and  the  subscript  i  indicates  the 
beginning  of  the  expansion  region.  Then  from  equations  (9)  and  (15): 
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This  approach  by  Johannesen  and  Meyer  can  be  simplified  by  the  use  of  a 
numerical  integration  technique  such  as  Runge-Kutta  to  solve  (<j>)  directly 
from  equation  (15), 

B.  The  Approaching  Flow  Solution. 


Beginning  with 

u  +  v  =0 
o  o 


(18) 


and  follow! rig  a  similar  procedure  the  following  equations  were  derived; 


u0  (4>)  "  wcosY  (19) 

v0  (4>)  =  wsinY  (20) 

ui  (4>)  =  Apcos  (2T  +  p-2)  +  Kj  (21) 

v1(<j>)  = -Aj  sin  (2T  +  82)  (22) 


where  4s~¥  =  constant 


and  A] .  Kq  and  82  are  constants. 


To  solve  for  the  initial  curvature  of  a  jet  boundary,  it  is  necessary  to 
consider  the  geometry  of  the  flow  at  the  lip,  see  Figure  2.  In  Figure  2, 


<t>fj  is  the  angl 
re spec 
as  R->0 


between  the  nozzle  wall  and  the  exit,  <j>L  and  <j>p  are  the 


re  angles  of  the  initial  and  final  Mach  lines  of 
Then 

the  expansion  region 

=  ^  -  PL  + 

(23) 

4>f  =  A  +  <t>L 

(24) 

4>j  =  p  +  4>f 

(25) 

A  =  6  +  UL  -  UF 

(26) 

where  and  bp  are  the  Mach  angles  bounding  the  expansion,  A  is  the  angle 
containing  the  expansion  and  5  is  the  angle  of  deflection  of  the  boundary 
streamline . 


The  boundary  conditions  along  the  wall  of  the  nozzle  are; 

uo  (4>N)  =  ~qL  (27) 

v0  (4>N)  -  0  (28) 

where  q  is  the  velocity  mangitude,  q  =-J  ujr  +  v^  in  the  uniform  flow  region. 
Then  by  equation  (20) 

Y  -  4>  -  4»N  ■  (29) 

Then  from  equations  (23)  through  (26) 
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Figure  2.  Flow  geometry  at  nozzle  exit. 


(30) 


u  (<i>T )  =  qL  cos  yL 

O  A-*  u 

vo(^L)  =  qL  Sin  Ul 

Combining  with  equations  (21)  and  (22)  gives 


(A)  =2  COS2]!!!,^)  -  sin  21-1^  C^)  +  cos  (<J>N> 


lqrL 


From  Figure  2  the  angle  of  the  streamline  is  given  by 


■Kl) 


0  .  4,  -  -  +  tan  xu; 

In  the  uniform  flow  region  8g>  -  0  and  the  boundary  conditions  at  t: 

5(j) 

at  R-0  is  v(0,<j)N)  v0(4>n>  =  0»  then  the  curvature  of  the  streamline  a 

origin  is  given  by, 

d9_  |  =  V1(<I>N) 

dR  R=0 


also  for  the  velocity ,  q 


q  dq 
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but  -U.  =  0- 

3<j) 
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and  equation  (32)  becomes 

U1  (^r)  =  qL  S±n  2U 


Tiie  geometry  also  gives 
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The  constant  is  determined  in  the  following  manners  From  equation 

(8) 

vq  =  u^  -  iL-  cos  [A(<j)  +  3o)]X 


then  dividing  equation  (31)  by  (30) 
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he  wall 
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(13)  and 
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tan  VL  =  YaJ&Jl 

u0 


then  equating  the  two  sets  of  relationships  for  the  velocity  components 

=  tan  yT 


vo(4>l) 

tan  = 

O 


from  which 


-♦l  +  iltan'1 


tan 


(39) 


and  equation  (4)  gives  a  relation  for  q: 

/9f\2 3 4 5 6  y-i  +  2  sin2  yL 

(  — )= - -  (40) 

\qL/  Y-l  +  2  sin2 

The  downstream  uniform  flow  field  imposes  the  boundary  conditions  that 
u0(j>j)  "  qF,  vo  ( 4> j )  =  0  so  that  (3p  =  4> j  and  the  final  velocity  in  the  expan¬ 
sion  region  is  qp.  Again  using  the  uniform  flow  equations  the  following  rela¬ 
tionship  for  the  initial  curvature  of  the  jet  boundary  is  derived: 


d9_i 

dR  R=0 


cot  yp 


dqj  — 
dR  r=q 


qF  sin  0_.  sin2  yF  + 
qp  sin  yF 


(41) 


The  Johannesen  and  Meyer  theory  can  be  used  in  plume  modeling  since  it 
allows  the  determination  of  the  initial  radius  of  curvature  of  the  plume  boun¬ 
dary.  The  next  section  will  illustrate  how  Korst  takes  advantage  of  this 
theory  to  match  plumes  from  two  geometrically  different  nozzles.  To  summarize 
the  procedure  for  obtaining  a  value  of  initial  curvature  of  the  jet  boundary 
the  following  steps  are  listed: 

(1)  The  flow  conditions  approaching  the  lip  are  known,  that  is  0^, 
dq  I  and  dSl 

dR  I  R=0  dR  1  R=0 . 

(2)  yL  and  up  (4>l)  are  calculated  from  the  known  flow  conditions  and 
equation  (37). 

(3)  qp  is  calculated  based  on  known  outside  pressure. 

(4)  yF  and  A  are  calculated  from  (40)and  (38)  through  (39) 
respectively. 

(5)  6  is  calculated  from  (30). 

(6)  The  integrals  of  equation  (16)  are  evaluated  between  the  limits 
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tanl-L  and  U  =  ^  +  AA. 

-U  Lj 

~T~ 


(7)  Equation  (16)  is  evaluated  for  uq  (4»p).  Or  equation  (15)  is 
integrated  numerically. 

(8)  Finally  equation  (41)  can  be  evaluated  to  obtain  the  initial  cur¬ 
vature  of  the  jet  boundary. 

III.  KORST' S  PLUME  MODELING  THEORY 


The  Johannesen  and  Meyer  theory  described  in  the  previous  section  provides 
a  method  for  determining  the  initial  radius  of  curvature  and  initial  slope  of 
the  jet  boundary  of  a  rocket  nozzle.  These  two  parameters  are  sufficient  to 
describe  the  initial  geometry  of  the  plume.  Korst-  expands  on  this  idea  by 
asserting  that  to  properly  model  a  prototype  nozzle  the  plume  geometry  of  the 
prototype  arid  model  must  be  the  same.  Korst  has  developed  a  numerical  intona¬ 
tion  technique  which  uses  the  theory  of  Johannesen  and  Meyer  to  first  deter¬ 
mine  the  radius  of  curvature  and  slope  for  a  specified  prototype  and  then  to 
determine  the  exit  flow  conditions  of  a  model  which  geometrically  matches  the 
prototype  plume.  It  is  generally  assumed  that  the  specific  heat  ratio  of  the 
model  and  prototype  are  not  the  same.  This  means  that  in  wind  tunnel  testing , 
air  or  some  other  inert  gas,  such  as  Freon,  can  be  used  rather  than  the  actual 
rocket  propellant. 


Consider  the  velocity  equations  of  Section  II  as  non-dimensionlized  by  the 
critical  speed  of  sound,  a*.  For  the  flow  approaching  the  nozzle  exit, 
equations  (19)  and  (20)  become  respectively? 


u 

0 

C4>l> 

a* 

=  cos 

V 

o 

(4>l> 

a* 

=  M£  sin  yL 

(42) 


(43) 


where  M*  is  the  critical  Mach  number.  Equation  (37)  which  gives  the  relation 
for  the  expansion  term  in  a  uniform  flow  region  may  be  written: 


U1^L') 

a* 


dM*, 
dR  ' 


R=0 


-  sin 


sm 


V 


(44) 


14 


;d0 

For  a  conical  nozzle,  0  =  constant,,  dR  =  0  so  that  equation  (44)  reduces  to 


4  sin  0,  cos  yT 
L  L 


2 

y-1-  i-xV 


(45) 


The  uniform  flow  approaching  the  nozzle  lip  is  also  the  initial  flow  for  the 
expansion  region.  The  constant  Cp  in  equation  (16)  can  be  evaluated  since  p  = 
Pl  so  that  the  integrals  vanish,  pl  can  be  determined  from  the  approaching 
flow,  4)^  can  be  evaluated  from  equation  (23)  so  that  (30  is  obtained  from 
equation  (39). 


Equation  (41)  can  also  be  re-written  in  terms  of  M*  as: 


2 

in0  sin  y  +  u  (<j>  )  1 

F  F  — —  M*  -I 

a*  F 

from  which  the  desired  result,  the  Initial  radius  of  curvature  of  the 
boundary,  rCj  can  obtained,  since: 

r  =  1  1 
c  i  d6 

dR  R=0 


d0 

dR 


R=0 


-  _1_  dM*  j 
*  dR 


-1 


“5 


coty  =  sin2y 


R=0 


(46) 


jet 


(47) 


The  assumption  of  a  conical  nozzle  also  defines  the  approaching  streamline 
angle  at  the  nozzle  wall  as  being  equal  to  the  conical  divergence  angle,  0^. 
Therefore,  the  Prandtl-Meyer  relations  can  be  used  to  find  the  final 
streamline  angle,  0p,  since: 

0F  =  0L  +U)  (V  (1V  (48) 

and 

4>f  =  ew  -  yp  +  a.  . 

F  F  2 

Geometric  modeling  of  the  plume  is  achieved  when  the  initial  slope  of  the 
jet  boundary  of  the  model  and  prototype  and  the  plume  radius  of  curvature  of 
the  model  and  prototype  are  the  same.  That  is  when: 

0FM  =  ®FP  (50> 

and 

rcM  =  rCP  .  <51> 
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It  is  assumed  that  for  the  prototype  all  the  specifying  parameters  are 
known.  For  the  model,  the  choice  of  a  propellant  gas  fixes  ym*  The  numerical 
procedure  used  by  Korst  is  essentially  as  outlined  at  the  end  of  Section  II  as 
simplified  by  the  conical  nozzle  assumption,  qp  =  0.  By  selecting  a  value  of 

dR 

Mlm  there  is  obtained  a  value  for  0^  since  by  equation  (50)  0p  ■  is  known  and 
equation  (48)  can  be  used  provided  Mp^  is  known.  By,  following  the  calculation 
procedure  a  value  for  r<^  is  obtained,  is  compared  with  a  value  of 


calculated  from  the  prototype  conditions,  the  value  of 


is  adjusted 


accordingly  and 


the  calculation  procedure  continues  until 


rcM  “  rcp* 


Before  the  iteration  process  can  begin,  Mp  ,  which  is  dependent  on  the 
external  flow,  must  be  determined  by  some  method.  Korst  refers  to  the  calcu¬ 
lation  of  this  parameter  as  the  closure  condition  for  the  wake.  There  are 
several  possibilities  for  the  choice  of  a  closure  condition  relating  to  the 
recompression  ratio  at  the  end  of  the  wake  or  to  the  conservation  of  mass  in 
the  wake;  however,  with  only  one  unresolved  parameter,  both  conditions  cannot 
be  satisfied  simultaneously.  Korst  chooses  to  match  the  invlscid  streamline 


deflection  pressure 


Yp  M2Fp 


rise  as  given  by 

mV 

M  M 


(52) 


for  a  weak  shock  recompression  or 


<Ym~1)  -  2YpM^p  -  (Yp-1) 

“yyT  yp  + 1 


(53) 


for  a  strong  shock  recompression. 

Since  Mp  is  dependent  on  the  jet  to  ambient  pressure  ratio,  the  model 
nozzle  is  designed  for  a  specific  Mach  number  and  altitude,  the  ’’design 
point".  This  procedure  calculates  a  value  for  the  model  nozzle  conical 
divergence  angle,  0p  ,  hut  does  not  determine  throat  radius  of  curvature  and 
nozzle  length.  The  method  utilized  to  calculate  these  parameters  is  discussed 
in  Section  IV. 


IV.  NOZZLE  S0ULTI0N  BY  METHOD  OF  DUTTON  AND  ADDY  WITH  METHOD  OF 
CHARACTERISTICS 

Korst* s  plume  modeling  theory  as  described  in  Section  III  provides  a 
method  of  determining  the  exit  Mach  number  and  conical  divergence  angle  for  a 
model  nozzle  which  has  the  same  initial  jet  boundary  radius  of  curvature  and 
slope  as  a  known  prototype  nozzle.  This  section  discusses  a  procedure  for 
determining  a  length  and  throat  radius  of  curvature  for  the  model  nozzle  which 
produces  the  exit  conditions.  The  procedure  utilizes  an  expansion  method 
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developed  by  Dutton  and  Addy^  to  solve  the  transonic  flow  in  the  region  of  the 
nozzel  throat  which  establishes  initial  conditions  for  a  method  of  charac¬ 
teristics  program.  The  method  of  characteristics  routine  solves  the  flow  in 
the  expanding  nozzle  until  the  Mach  number  at  the  wall  boundary  is  equal  to 
the  final  Mach  number,  Mp^,  determined  by  the  modeling  program. 

Because  the  equations  for  the  method  of  characteristics  are  dependent  upon 
the  ralationship  ^Mz-1,  other  methods  must  be  used  to  solve  the  transonic 
flow  region  in  convergent-divergent  nozzles.  Hall^  developed  a  widely  used 
expansion  technique  which  is  a  small-perturbation  method  from  the  one¬ 
dimensional  flow  solution.  Hall's  solution  in  cylindrical  coordinates  is 
given  by 

u  =  1  +  ui  (r,  z)  e  +  U2  (r,  z)  e  2  +  U3  (r,  z) e  3  +  .. .  (54) 


vq  (r,  z)  £  +  V2  (r,  z)  £  2  +  V3  (r,  z)  £ 3  +  . . . 

mi 


(55) 


where  u  and  v  are  axial  and  radial  velocity  components  normalized  by  the 
sonic  velocity,  a*.  The  expansion  variable  £  is  defined  to  be  1/R.  R  is  the 
normalized  throat  radius  of  curvature,  R=Rw/r*,  where  Rw  is  the  actual  throat 
radius  of  curvature  and  r*  is  the  nozzle  throat  radius.  The  transformed  nor¬ 
malized  axial  coordinate,  z,  and  the  normalized  radial  coordinate,  r,  are 
defined  by  z  =/2R  \  /2  j*_and  r  =  y  . 

\y+l/  r*  j  r* 

Hall's  equations,  (54)  and  (55)  are  well-behaved  provided  R>1.5.  Kliegel 
and  Levine3  had  seemingly  overcome  this  limitation  by  using  a  series  expansion 
where  £  =  1/(R+1)  for  axisymmetric  nozzles.  Kliegel  and  Levine  developed 
their  series  solution  in  torodial  coordinates  such  that  the  axis  is  repre¬ 
sented  by  a  line  t|  =  0  and  the  nozzle  wall  is  represented  by  a  line  n  =  %• 

The  results  were  then  transformed  back  to  cylindrical  coordinates. 
Unfortunately,  the  series  do  not  satisfy  the  governing  differential  equations 
of  motion  in  cylindrical  coordinates.  This  paper  by  Kliegel  and  Levine3  is 
still  important  in  that  the  authors  re-derived  Hall's  original  equations  and 
corrected  errors  in  the  third  order  terms. 

Dutton  and  Addy^  formulated  the  problem  in  a  similar  manner  to  Hall  except 
that  the  expansion  veriable,  £  =  l/(R+ri)  was  used.  Details  of  the  series 
derivation  can  be  found  in  references  4  and  6  and  will  not.be  presented  here. 
The  basic  approach,  however,  is  to  substitute  equations  (54)  and  (55)  into  the 
governing  equations  and  into  the  boundary  conditions.  The  governing  equations 
are  taken  to  be  the  gas  dynamic  equation  and  the  irrotationality  condition. 

The  boundary  conditions  are  that  both  the  axis  of  symmetry  and  the  nozzle  wall 
are  streamlines.  After  substituting,  coefficients  of  like  powers  of  E  are 
gathered  to  formulate  the  various  orders  of  the  expansion.  The  equations  for¬ 
mulated  in  this  manner  are  then  solved  by  assuming  solution  forms  suggested  by 
the  boundary  conditions. 
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Dutton  and  Addy  suggests  using  third  order,  p  =  1  solutions  for  most 
applications.  This  is  the  approach  taken  in  the  program  presented  in  this 
report.  The.  equations  for  the  third  order  velocity  terms,  the  non-dimensional 
Mach  number  and  the  streamline  deflection  angles  are  given  as  follows; 


u (r,z)  -  1/2  f2  “  Vk  +  2 
u2  (  r ,  z  )  =  2y-f9r4  -  4y+15-12ri.  r2  +  10y+57-72T)  + 


24 


24 


288 


(56) 

2y-3  ^  z2  (57) 
6' 


u3  (r  »z)  =  556y2+1737y+3069r6  -  388y2  +  (1161-384t))  y  +  (1881-1728n)  r4 
10,36.8  2304 

+  304y2  1-  (831-<576ri)  Y  4-  (1242-2160ti+864r)2 )  r2 
1728 

9 

^60ri)  y  +  14, 211-32. 832r)>  -f  20,736p2) 

82,944. 

+  |j52yZt5lY+327  r4  -  52Y24-(75y2+(279-288n)r2+  ?2j_2+18_0y  +(639-1080r]+432n2)  |  z 
384  r  192  1152  J 


f  f»2y^3  r2  +  ( 1 3-1 6'f) ) y- ( 2 7-2 4 rj )~j  r2  +  [4y2-i57Y+27  ]  r3 

L  R  . 48  —  “J  ! 14 5— 


&  48 

vi  (v,  •.)  =  l/4r3  -  l/4r  +  rz 


96 


288 


V',*  (r9z)  =  y+3  rD  -  20y-f63--36r)  rJ  +  28y+93-108'n  r  +  ^2y+9  -  4y+15-12ri  r|  z 


L.Q  —3 


(58) 

(59) 
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+  rz' 


(60) 


v3  (r  sz)  =  6836y2+23.031vf30,627  r7  3380y2-Hil.391-3840n)y-K15,291-ll,520n)  r5 
82,944  13,824  '  . 

+  3424y2  +  (UA271-7200ri)Y+(15,228-g2,680'n+6480ri2)r3 

13,824 

-  710Qy2  +  (22,311-2Q,160ti)y  +  (30, 249-66, 960p+38, 880ri2 )  r 


82,944 

+r556Y2+1737Y+3069  r5  388Y2+(1161-384q)y-K1881--1728n)r3 
L  “  1728  576 

-f304y2+(831-576ri)Y+(1242-2160q  +  864ri2)  rlz 
864  J 


■f52Y2+51y+327  r3  -  52y2-!-75Y+(279-288q)^  z  -t  f-  7y-3  rl 
L  192  192  J  L  12  J 

M*  (r  >  z)  =1  +  u^e  +  u^2  u3  -I-  “T~  v2  ^  e3  +  ... 


(61) 


(62) 
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~Y+l 

% 

0  (r , z)  = 

2  e 

£  +  <V2“U1V1)  e?+  (V3“U1V2_U2V1+U1V1>  6  3 


+  .  .  . 


(63) 


This  result  is  used  by  assuming  that  u»v  tnd  the  streamline  deflection 
angle,  0,  is  small  then: 


M*  -  u  (64) 

0  «  v  (65) 

since  v  is  non-dimensionalized  by  a*  and  in  the  region  of  the  throat  a*=M*. 

0  is  given  by  tan  0  =  v/u  and  for  small  0i  0  =  =  v/u.  The  calculation  proce¬ 
dure  begins  with  the  user  choosing  a  value  for  R.  The  position  on  the  axis 
where  conical  wall  joins  the  circular  throat  can  be  calculated  from: 

Xc  -  Rsin0c  (66) 

where  ©c  is  the  conical  divergence  angle  of  the  nozzle  and  Xc  is  non- 
dimensionalized  by  r*.  The  distance  from  the  axis  to  the  wall,  Y,  at  any  axis 
position,  X,  is  given  by: 

Y  =  1  +  R  (1  -Jl  -  x  2 

^  R  (67) 

Figure  3  illustrates  this  relationship.  The  velocity  components  as  given  by 
(54)  and  (55)  are  evaluated  at  increasing  X  until  the  axis  velocity  is 
slightly  supersonic.  At  this  point,  values  of  u  and  v,  hence  M*  and  0,  are 
calculated  from  the  axis  to  the  nozzle  wall. 

The  values  of  M*  and  0  along  with  the  X  and  Y  positions  become  the  initial 
values  for  a  method  of  characteristics  routine  which  solves  the  flow  in  the 
nozzle  until  the  local  Mach  number  at  the  nozzle  wall  matches  the  model  exit 
Mach  number  determined  from  the  modeling  theory.  Then  the  X  position  at  this 
point  is  the  model  nozzle  length. 

V.  PLUME  MODELING  EXPERIMENTS 

The  Aeronautical  Research  Institute  of  Sweden  (FFA)  has  developed  a  hot 
gas  system  for  the  FFA  .5  x  .5  meter  supersonic  wind  tunnel  for  the  purpose  of 
studying  plume  effects.  The  efforts  of  FFA  have  been  closely  coordinated  with 
members  of  the  Gas  Dynamics  Laboratory  at  the  University  of  Illinois  at 
Urbana-Champaign. 

Experiments  were  conducted  by  FFA  with  the  expressed  purpose  of  critically 
evaluating  the  merits  and  limitations  of  plume  modeling  techniques.  Some  of 
this  effort  has  been  documented  by  Ifyberg  et.  al.^.  For  the  plume  modeling 
experiments  reported  in^,  FFA  designed  two  air  nozzles  to  be  used  as  the  pro¬ 
totype  nozzles  and  two  model  Freon  nozzles,  one  for  weak  shock  modeling  and 
one  for  strong  shock  modeling.  The  geometrical  matching  of  plume  shapes  was 
verified  by  Schlieren  photographs.  Pressure  measurements  were  made  to  verify 
the  matching  in  the  base  region. 
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The  matching  of  the  plume  shapes  as  presented  in  reference  7  is  shown  in 
Figures  4a  and  4b.  Figure  5,  also  from  reference  7,  demonstrates  the  base 
pressure  matching  where  P^  is  the  base  pressure,  P^  is  the  exit  pressure  and 
P  is  the  pressure  at  the  nozzle  lip.  From  these  experiments,  FFA  concluded 
tnht  the  base  pressure  agreement  was  satisfactory  at  both  the  design  point  and 
for  a  wide  range  of  off-design  conditions.  The  overall  conclusion  of  FFA  was: 
"The  Freon  plumes  shapes  have  been  found  to  be  in  close  agreement  with  those 
of  the  corresponding  air  test  supporting  the  suggested  modeling  methodology 
and  design  precedures . Testing  by  the  Army,  which  will  be  documented  at  a 
later  date  tends  to  confirm  this  observation. 
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©  AIR  PROTOTYPE 

ID  FREON  WEAK  SHOCK  MODELING 

A  FREON  STRONG  SHOCK  MODELING 


APPENDIX  I 


SAMPLE  CASE 

Tables  I  and  II  are  the  ordinary  output  from  the  model  design  program. 
Much  more  detailed  nozzle  flow  field  information  and  de- bugging  information 
can  be  obtained  by  activating  logical  unit  6.  There  are  some  Perkin-Elmer 
Interdata  8/32  dependent  statements  in  the  program.  For  instance,  subroutine 
PRTOPT  which  assigns  the  logical  units  for  output  and  function  IGC  which  is 
part  of  the  free  formatted  read  subroutine  would  have  to  be  re-written  for 
use  on  any  other  computer.  However,  the  bulk  of  the  program  is  written  in 
standard  Fortran  IV. 

The  program  is  currently  written  to  accept  inputs  interactively  from  a 
terminal.  The  prototype  inputs  for  the  sample  case  were:  specific  heat 
ratio,  nozzle  angle,  exit  Mach  number  and  jet  surface  Mach  number.  Model 
inputs  were:  specific  heat  ratio  and  throat  radius  of  curvature.  All  the 
other  parameters  shown  in  Table  I.  were  calculated  by  the  program.  The  weak 
shock  modeling  condition  was  chosen  for  the  sample  case. 

Table  II.  tabulates  the  calculated  plume  shapes  for  the  model  and 
prototype.  The  prototype  plume  shape  was  calculated  using  the 
Johannesen-Meyer  theory  and  the  model  plume  shape  was  calculated  using  a 
method  of  characteristics  routine  which  is  started  with  exit  conditions 
calculated  from  a  separate  nozzle  method  of  characteristics  routine.  It  is 
seen  from  Figure  6.  that  there  is  excellent  agreement  between  the  two  plume 
shapes  for  approximately  one  body  diameter. 
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MODEL  NOZZLE  DESIGN 


PROTOTYPE 

SPECIFIC  HEAT  RATIO. ............................  1.235 

NOZZLE  ANGLE ....................................  1 0 . 360 

EXIT  MACH  NUMBER. . ..............................  2.530 

JET  SURFACE  MACH  NUMBER .........................  3-350 

INITIAL  SLOPE  OF  JET  PLUME. .....................  32.304 

INITIAL  RADIUS  OF  CURVATURE  OF  JET  PLUME........  5.426 

PRESSURE  RATIO. ................................ .  0.012 


MODEL 

SPECIFIC  HEAT  RATIO.  ............................  1.400 

THROAT  RADIUS  OF  CURVATURE. . ....................  3.000 

BEGINNING  AXIAL  LOCATION  OF  CONICAL  SECTION.....  0.207 

BEGINNING  RADIAL  LOCATION  OF  CONICAL  SECTION....  1.00? 

NOZZLE  ANGLE,  DEG...............................  3.964 

NOZZLE  LENGTH...................................  2.952 

NOZZLE  EXIT  RADIUS. .............................  1.197 

EXIT  MACH  NUMBER ................................  1.761 

JET  SURFACE  MACH  NUMBER. ........................  2.90? 


TABLE  I.  MODEL  AND  PROTOTYPE  CHARACTERISTICS. 
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PLUME  SHAPE  CALCULATED  FROM  JOHANNESEN-MEYER  THEORY 


X 

R 

THETA 

0.0000 

1 .0000 

32.3042 

0.1000 

1.0617 

31 .0633 

0.2000 

1.1205 

29.8384 

0.3000 

1.1765 

28.6283 

0.4000 

1 ,2297 

27.432 

0.5000 

1 .2803 

26.2486 

0.6000 

1.3283 

25.0771 

0.7000 

1.3739 

23.9167 

0.8000 

1 .4170 

22.7667 

0.9000 

1.4578 

21.6262 

1 .0000 

1.4964 

20.4947 

1.1000 

1 .5326 

19-3714 

1.2000' 

1.5667 

18.2559 

1.3000 

1.5986 

17.1474 

1.4000 

1.6284 

16.0456 

1.5000 

1 .6561 

14.9498 

1.6000 

1.6818 

13-8596 

1.7000 

1,7055 

12.7744 

1.8000 

1.7272 

11.6939 

1.9000 

1 .7469 

10.6176 

2.0000 

1 .7647 

9.5451 

2.1000 

1 .7805 

8.4760 

2.2000 

1 .7945 

7.4098 

2.3000 

1.8065 

6.3462 

2.4000 

1 .8167 

5.2848 

2.5000 

1 .8250 

4.2252 

2.6000 

1.8315 

3.1670 

2.7000 

1.8361 

2.1099 

2.8000 

1 .8389 

1.0536 

2.9000 

1 .8398 

-0.0024 

3.0000 

1 .8388 

-1.0585 

MODEL  PLUME  SHAPE  CALCULATED  BY  METHOD  OF  CHARACTERISTICS 


X 

R 

THETA 

0.0000 

1 .0000 

32.3036 

0.4153 

1.2387 

27.5741 

0.6224 

1.3430 

25.7525 

0.8406 

1.4446 

24.0213 

1.0653 

1 .5411 

22.4039 

1.0654 

1 .5412 

22.4048 

1.3355 

1 .6478 

20.5551 

2.6375 

2.0424 

13.1913 

2.6490 

2.0451 

13.0891 

TABLE  II.  MODEL  AND  PROTOTYPE  PLUME  SHAPES 
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APPENDIX  II. 

MODEL  NOZZLE  DESIGN  PROGRAM  LISTING 


C«  «  «4«»  *»«*»»*»«#»»» 

C 

C  MODEL  DESIGN  PROGRAM 

C 

Qi't»sft«sa!4{sa!i{!j(»6i}eiJsa»sa8s»#tSM#»a9a»jH>#«s«8a  '4;  ssa»#ess«aa{ia 

C 

REAL  ML  P , MFP , MFM , MLM 
INTEGER  Z ,  ZO ,  S ,  P ,  Q ,  A8 

DIMENSION  XA(3,45) ,R(3,45) ,XM( 3 ,45) ,T(3 , 45) , 

»HOL(20) ,FL0(20) ,B(45) ,0(45) 

COMMON/CKLEV/B ,C 

COMMON/SHAPE/X7 , CO , R7 , TO , RO ,GAMM,MLM 
COMMON/M ARRAY/XA , XM , R , T 
COMMON/INDICE/I, N,K, IS, IP 
COMMON/I NJET/MEXIT , MFM , KO , NO 
COMMON/VAR/VAR(17) 

C 

5  CALL  PRTOPT 
C 


C 

CALL  MODEL( MFM ,MEXIT , THTALM) 

C 

CS»S89#8881i9«»«»***»«89«###»S#8»«88»»«»9*«»»a888*8a»»8»9»S8»»»«9«a# 

Q*«*«»*k**a»»»**»«»«*»X»*»*tt*»»»»t«***»**«tt**ft»«*«**#*«**»ft»*»*tttt»* 

C 

CALL  NOZZLE(MEXIT, THTALM) 

C 


Caaft»8a#fta»»»»ait«tt8*»4>«»««8tttf  »*»«*««»»»«««•& «*«»»«*«»» 

C 

CALL  PLUME 
C 

QaaaaaaaBftasaaaaaaaBaBBBttttaBaaBsagaaaaaaaaattaaaaflaaaaaaaaaaaaffaaaBB 

C 

CALL  PLUL1S 
C 

C*«»*»tt**«***ii»**«***»***«#»»*»»*»«fttt»««tf«**»*«**«*tf**»*»»**K«*»**« 
£****» *****»»*****»»»»**«***»*«*«*««* *««**»*»»«{***«**»*#«»*»******» 
C 

C  END  OPTIONS 
C 

WRITE (1 ,130) 

130  FORMAT ( /T20' ENTER : ' 

*T27'(1)  RUN  ANOTHER  CASE'/ 

*T27' (2)  STOP' ) 

CALL  FFREAD(HOL, FLO,LH,LF) 

IF( FLO( 1 ) . EQ. 1 . )  GO  TO  5 
C 

STOP 

END 
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SUBROUTINE  MODEL ( MFM , MEXIT , THTALM) 

C 

C  CONTROLS  THE  MODELING  SECTION  OF  THE  PROGRAM 
C 

REAL  MLP , MFP , MFM , MLM 
INTEGER  Z ,  ZO ,  S ,  P ,  Q ,  A8 

DIMENSION  XA(3r45)',R(3.45)  ,XM(  3  » 45) ,  T(3 , 45)  ,HOL(20)  ,FL0(20)  , 
*B(45) ,C(45) 

COMMON/ CKLEV /  B ,  C 

COMMON/SHAPE/X 7 » CO , R7 » TO , RO , G AMM , MLM 
COMMON/MARRAY/XA ,XM ,R,T 
C 

COMMON/RADS/PI, RAD 
C 

PI=3* 141592654 
RAD- 5 7 • 3 
1  IERRsO 

C 

C  SPECIFY  TYPE  OF  INPUT 
C 

WRITE ( 1 j 10) 

10  FORMAT (/T31 ENTER : * 

S/T5 ' 1  PROTOTYPE  INITIAL  JET  SLOPE  AND  RADIUS  OF  CURVATURE8 
ss  KNOWN' 

*/T5'2  PROTOTYPE  NOZZLE  GEOMETRY  AND  EXIT  CONDITIONS  KNOWN*) 

C 

CALL  FFREAD(HQL,FLO,LH,LF) 

INOPT=FLO( 1 ) 

C 

GO  TO  (20,30) , INOPT 
C 

C  PLUME  SHAPE  TO  BE  SPECIFIED 
C 

20  CALL  PLUMIN( THTAFP , RADP) 

GO  TO  40 
C 

C  PLUME  FLOW  CHARACTERISTICS  TO  BE  SPECIFIED 
C 

30  CALL  FLOWIN(G AMP ,  MLP ,  THTALP ,  MFP ) 

C 

C  SPECIFY  MODEL  FLOW  CONDITIONS 
C 

40  CALL  MODIN (GAMM, MFM, GAMP, MFP, IERR) 

IF(IERR.EQ. 1 )G0  TO  1 
C 

C  CALCULATE  PLUME  SHAPE  IF  PROTOTYPE  FLOW  SPECIFIED 
C  (IN0PT=2)  AND  PRINT  RESULTS 
C 

IF(INOPT.EQ. 1 )G0  TO  50 

CALL  CALCPL(GAMP , MLP , THTALP , MFP , THTAFP , RADP ) 

CALL  LISTPR( MFP, GAMP) 

C 

C  BEGIN  ITERATION 
C 

50  CALL  HERAT ( GAMM, THTAFP, THTALM, RADP, MFP, MFM, MLP, MLM) 

C 


30 


RETURN 

END 


SUBROUTINE  NOZZLE (MEXIT, THTALM) 

C 

C  CONTROLS  THE  SECTION  OF  PROGRAM  WHICH  CALLS  NOZZLE  SOLUTION 
C  AND  USES  METHOD  OF  CHARACTERISTICS  TO  CALCULATE  THE  NOZZLE 
C  LENGTH . 

C 

REAL  ML P  8  MFP , MFM , MLM 
INTEGER  Z,Z0,S,P,Q,A8 

DIMENSION  XA(3845) ,R(3,45) ,XM(3s45) ,T(3,45) 

»,B(45) ,C(45) » HQL ( 2  Q ) , FLQ ( 2  Q ) SQA(45) ,PA(45) 

CQMMON/CKLEV/B ,C 

COMMON/SH APE/X7 , CO » R7 , TO , RO , GO , XMO 
COMMON/MARRAY/XA ,XM , R , T 
COMMON/VAR/VAR(17) 

TO=THTALM®57. 3 

Ks9 

C 

MEXIT= 1 
VAR(8)=G0 
VAR( 1 5) =XM0 
VAR( 1 2) sTO 
TOsTO/57- 3 

74  WRITE ( 1 ,75) 

75  FORMAT); /T20 8 ENTER  THROAT  RADIUS  OF  CURVATURE’) 

CALL  FFREAD(HOL,FLO,LH,LF) 

CO=FLO(1) 

VAR(9 )~C0 
C 

WRITE (1 ,138) 

138  FORMAT(/T20» ENTER:  1 -NOZZLE  SOLUTION  WITH  KLIEGEL-LEVINE * 
1,/T27'2-NOZZLE  SOLUTION  WITH  ADDY-DUTTQN’ ) 

CALL  FFREAD(HOL,FLO,LH,LF) 

ISUB=FLO( 1 ) 

X?-COsSIN(TO) 

R7  =  1 .+C0*d .-SQRT(1 .-(X7/CO)**2) ) 

VAR(10)=X7 
VAR( 1 1 ) =R7 
C 

WRITE(6,300)K 

300  FORMAT ( /IX /NUMBER  OF  INITIAL  POINTS  ='I3) 

C 

C  CALL  TO  NOZZLE  SUBROUTINES  TO  SOLVE  THE  FLOW  FIELD 
C  CLOSE  TO  THROAT  CONTINUES  UK! EL  X  POSITION  REACHED  WHERE 
C  MACH  GREATER  THAN  1.025. 

C 

C 

DO  90  X=0.0,X7,(X7/10.)-. 000001 
Y0= 1 . +C0* ( 1 ,-SQRT( 1 . - (X/CO) a#2) )  . 

C 

J=1 

DO  100  Y=0,Y0,(Y0/FLOAT(K-1)) 

IF(ISUB.EQ. 1 )CALL  KLGLEV(X , Y , J) 

IF(ISUB.EQ.2)C ALL  ADDUT(X,Y,J) 

PA( J )sX 
QA( J ) sY 

WRITE (6 , 400)X , Y 
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400  FORMAT(1X, 'X  =  ' , F20.8, '  Y  =  ' ,F20.8) 

J  =  1 

100  CONTINUE 

IF(B(1).GT.1.025)G0  TO  32 
90  CONTINUE 

WRITE(1 ,401 ) 

401  FORMAT ( /T5 ' **  SUPERSONIC  FLOW  NOT  ACHIEVED  INCREASE  THROAT’ 
a'  RADIUS  OF  CURVATURE  **') 

WRITE(1 ,402)B(1) 

402  FORMAT(/T5 'HIGHEST  MACH  BEFORE  CONICAL  SECTION  =  ’F8.4) 

GO  TO  74 

C 

C  SETS  UP  INITIAL  CONDITIONS  FOR  THE  METHOD  OF  CHARACTERISTICS 
C  PROGRAM  WHERE  XA  ARE  THE  X  LOCATIONS  (ALL  ARE  EQUAL  INITIALLY) 

C  R  ARE  THE  RADIAL  LOCATIONS 

C  XM  ARE  THE  DIMENSIONLESS  SPEED  MACH  NUMBERS  M* 

C  T  ARE  THE  THETAS  (ANGLE  BETWEEN  THE 

C  X  CO-ORDINATE  AND  THE  VELOCITY  VECTOR. 

C 

32  DO  110  J=K  ,1,-1 

XA(1 ,K+1-J)=PA(J) 

R( 1 ,K+1-J)=QA(J) 

XM( 1 ,K+1-J)=B(J) 

T(1 ,K+1-J)=C(J) 

110  CONTINUE 
C 

C  INITIALIZE  POINTERS  USED  IN  THE  METHOD  OF  CHARACTERISTICS 
C 

1=1 
S=  1 
P=1 
A8=2 
Q=0 
N  =  0 
Z0=0 
KO=K 
IFINI=0 
C 

DO  304  J=1,K 

WRITE(6,303)J»XA(rfJ)fR(I,J) ,XM(I,J),T(I,J) 

303  FORMAT(1X,I3,2X,4F8.3) 

304  CONTINUE 
C 

WRITE(6 , 305) 

305  F0RMAT(/,5X, ' FLOWFIELD  COMPUTED  ALL  POINTS  PRINTED' ) 

C 

C  BEGIN  METHOD  OF  CHARACTERISTICS  LOOP 
C 

1260  CALL  SE'TPTS(K ,  S ,  I) 

IF(I.NE.2)G0  TO  1320 
C 

CALL  WALL(I) 

CALL  CNTLIN(K , S, I) 

CALL  YINCR(P,MEXIT,I,K,S) 

1320  CALL  SETEXT(P , A8 , I ,Q , K , S, N , IFINI ) 

IF(IFINI.EQ. 1)G0  TO  140 
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ITEST-(I/2)B2 


IF(ITEST.NE.'I)A8=2 
CALL  CROSSd jKjS, A8) 


1=1+1 

S=S+1 

IF(I.NE.3)GO  TO  1260 
CALL  RESET (K, S) 

1-1 

GO  TO  1260 

RETURN 

END 


SUBROUTINE  PLUME 


C 

C  CONTROLS  THE  SECTION  OF  THE  PROGRAM  WHICH  CALCULATES 
C  THE  MODEL  PLUME  SHAPE  USING  METHOD  OF  CHARACTERISTICS 
C 

REAL  MLP , MFP , MFM , MLM 
INTEGER  Z,ZO,S,P,Q, A8 

DIMENSION  XA(3,'45),'R(3,45),XM(3,45)',T(3,45) 

*,B(45) ,C(45) 

COMMON/ CKLEV/B ,C 

COMMON/SHAPE/X7 , CO , R7 , TO , RO , GAMM , MLM 
COMMON/MARRAY/XA , XM ,  R , T 
COMMON/INDICE/I , N , K , IS , IP 
COMMON/IN JET/MEXIT , MFM , KO , NO 
C0MM0N/VAR/VARO7) 

C 

C  CALCULATE  THE  ACTUAL  PLUME  SHAPE  USING  METHOD  OF  CHARACTERISTICS 
C 

N0=9 

KO=NO 

K=KO 

N=N0 

CALL  JINIT 
1=1 

270  CALL  SETPT2 
CALL  WRITRW 
1=1+1 
IS=IS+1 

IF(I.GT.(N-1))G0  TO  1210 
GO  TO  270 
C 

1210  J=1 

1215  CALL  PLUMPT 
CALL  SETPT3 
N=N+1 
IP=IP+1 
C 

C  CHECK  FOR  END  CONDITION 
C 

IF(N.GT.(N0+K0-2) )RETURN 
K=K-1 

GO  TO  1215 
C 

RETURN 

END 
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SUBROUTINE  ADDUT(X,Y,J) 

REAL  N,N2 

DIMENSION  B(45)SC(45) 

COMMON/ CKLEV/B , C 

COMMON/SHAPE/X7 ,  CO ,  R7 ,  TO ,  RO ,  G ,  XMO 
C 

Nsl.  - 

E=1 ./ (CQ+N) 

Z  =  X/SQRT(((G+1.)/2.)*E) 

Y  2=Y#Y 
Y3-Y2*Y 
I4sY3#Y 
Y5=Y4*Y 
Y6sY5aY 
Y?=Y6SY 
C 

.  E2=E*E 
E3-E2«S 
C 

Z2-Z»Z 

Z3-Z2«Z 

C 

G2-G*G 

C 

N2=NaN 

C 

U1=.5eY2-.25+Z 
V1  =  (  . 25®Y2~.  25-s-Z)  SY 

U2s( (2.BG+9. ) /24.)*Y4  -  ( (4 .aG+15.-12.8N) /24. ) ®Y2 
A  +  (10. *G+57 . -72 . 8N ) /288 .  +  (Y2+(4.*N-5. ) /8.)*Z 
B  -  ((2.*G-3.)/6.)*Z2 
C 

V2=((G+3.)/9  J*Y5  -  ((20.sG4-63.”36.sN)/96.)«Y3 
A  +  ( (28. sG+93 •- 108 .sN)/288.)aY 

B  +(((2.*G+9.)/6.)*Y3  -  ((4.*G+15.-12.*N)/12.)*Y)*Z+Y*Z2 
C 

U3-  ( (556  .  'SG2+1737  .aG+3069 . )  /  1 0368. )  SY6 
A  -( (388.sG2+(1 161 . -384 . 8N) *G+( 1 88 1 . - 1 728. ) )/2304 . ) *Y4 

B+((304.aG2+(831 .-576. aN)aG+( 1242. -2160. «N+864.*N2))/1 728.) »Y2 
C- (2708. sG2+( 7839. -5760.aN)aG+( 14211 .-32832. aN+20736.sN2) )/82944. 
D+(  (52.aG2-«-51  .aG+327. )  /384. )  aY4*Z 
E-( (52.aG2+75 .aG+(279 .-288 . 9N ) )/1 92. ) aY2*Z 
F+(  (92.&G24-l80,»G+(639.-1080  .»N4-432.‘sN2)  )/1 152. )*Z 
G-( (7.aG-3.)/8.) »Y2*Z2 
H+(((13.-l6.*N)*G-(27.-24.»N) )/43.)aZ2 
I+( (4 .*G2-57 -*G+27. ) /1 44. ) *Z3 
C 

V3=(  (6836.9G2+23031  .aG-f30627.)/82944.)»Y7 
A-( (3380. aG2+(1 1391 .-3840.«N) aG+( 15291 .-11520.2N))/1 3824.) »Y5 
B+((3424.*G2+(1 1271 .-7200. aN) »G+( 1 5228 . -22680. 3N+6480.aN2) 

1  ) /I  3824 , ) ®Y3 

C-(  (7  100  .IJG2+(2231 1 .  -20160 .  3N)  *G+(30249 .  -6696O .aN+38880 .aN2) 

1  ) Z82944 . ) aY 

D+((556.aG2+1737.sG+3069.)/1728.)«Y5j»Z 

E-((388.*G2+0 161 .-384. sN)3G+( 1881 .-1728.*N) )/576.)»Y3»Z 

F+( (304.aG2+(831 .-576. aN)»G+( 1242. -2160. aN+864 .aN2) )/864 . ) aYaZ 
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o  o 


G+( (52.»G2+51 .aG+327.)/192.)*Y3*Z2 
H-( (52.*G2+75.sG+(279.-288.*N) )/192.) 4Y*Z2 
I-((7.*G-3.)/12.)*Y8Z3 
C 

B(J)  =  1  .+U1«E+U2*E2+((U3+(G+1  .)/.25)*V1*'/1)*E3 
C 

C(J)=SQRT((G+1 .)/2.*E)*(V1*E+(V2-U1*V1)*E2 
A  +( V3~U 1*V2-U2*V1+U 1*U 1*V T)*E3) 


RETURN 

END 
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o  o  o  o  o  o 


SUBROUTINE  CALCPL(GAMMA , MACHL , THTAL, MACHF . THTAF , RADIUS ) 
REAL  LAMDA , MUL , MACHL , MACHF ,MUF 
COMMON/RADS/P I , RAD 
C 

C  THIS  ROUTINE  CALCULATES  THE  INITIAL  SLOPE  OF  THE  PLUME 
C  AND  THE  INITIAL  RADIUS  OF  CURVATURE  OF  THE  PLUME. 

C  THE  FLOW  OUTSIDE  OF  THE  CENTERED  WAVE  REGION  IS  CALCULATED 
C  BY  ESTABLISHING  BOUNDARY  CONDITIONS  FROM  MATCHING 
C  SOLUTIONS  AT  THE  LOWER  (PHIL)  AND  FINAL  (PHIF)  MACH  LINES 
OF  THE  CENTERED  FAN  REGION. 


LAMDA=SQRT( (GAMMA- 1 . ) /(GAMMA+1 . ) ) 

MUL=XMANG 1(1. /MACHL ) 

PHIL-THTAL-MUL+FI/2 . 

E1sATAN(LAMDA/TAN(MUL) ) 

BO  =E 1 /L  AMD A- PHIL 

XI  AND  X2  ARE  THE  M»  MACH  VALUES  X1=M*L  AND  X2=M*F 

X 1 sXMSQR ( MACHL , GAMMA ) 

X2-XMSQR ( MACHF , GAMMA ) 

C 

C  CALCULATES  THE  PRANDTL-MAYER  FUNCTIONS  FOR  M*L  AND  M«F 
C  XL1=W(M*L)  AND  XL2=W(M*F) .  THEN  XL2-XL1  IS  THE 
C  DELTA  TURNING  ANGLE  DUE  TO  THE  LIP. 

C 

CALL  TANGLE (XI, XL  1, LAMDA) 

CALL  TANGLE ( X 2, XL2, LAMDA) 

C 

C  THE  FOLLOWING  EQUATIONS  CALCULATE  THE  APPROACHING  FLOW 
C  IN  REGION  B  (OUTSIDE  OF  THE  CENTERED  WAVE) 

C 

C0=SIN(THTAL)^(4 (CQS(MUL) )«#2/ 

S(MACHL'SMACHL-1 . ) - (SIN(MUL) ) a82 ) 

C 

C2=-2 .®LAMDA*SQRT(X 1 )/( COS(£1 )S8( ( (3 .®GAMMA)- 1 . ) /(2 .* (GAMMA- 1 . ) ) ) 
1sSIN(E1 )S9 .5) 

C 

C1=C0gC2 

E2=E 1 +LAMDA# (XL2-XL 1 -XMANG 1 ( 1 . /MACHF )+XMANG1 ( 1 ./MACHL) ) 

C 

C  THE  FOLLOWING  SECTION  INTEGRATES  THE  INTEGRALS  OCCURING 
C  IN  THE  CENTERED  'WAVE  EQUATIONS. 

C 

D1=(E2-E1)/10. 

C 

S1=0.0 
S2=0 . 0 
.  S3-0.0 

S4=0 .0 
C 

DO  10  XN=1 .  ,10. 

E3=E1+D1*(XN-.5) 

S1=C0S(E3/LAMDA)8SIN(E3)®s(-.5)*C0S(E3)8* 

*(-1./(2.*LAMDA*LAMDA))*D1+S1 
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c 

S2=SIN(E3/LAMDA)*SIN(E3)ss(™.5)*COS(E3)*« 

*(-1 ./(2.*LAMDA*LAMDA))*D1+S2 
C 

S3=SIN(E3/LAMDA)#SIN(E3)*»(-1 .5)*COS(E3) 

*( (GAMMA-3. ) /(2. * (GAMMA- 1. ) ))«D1+S3 
C 

S4=COS(E3/LAMDA)*SIN(E3)"(-1 .5)aGOS(E3)#a 
*( (GAMMA-3. ) /(2 .* (GAMMA-1 . ) ) )®D1+S4 
C 

10  CONTINUE 
C 

U0=-((C0S(E2)8*((3.*GAMMA-1 .) /(2 .* (GAMMA 
*-1.) )))*SIN(E2)**.5)/(2.*LAMDA) 

C 

c 

C  U1  IS  THE  VELOCITY  AT  THE  UPPER  LIMIT  OF  THE  CENTERED  WAVE 
C  REGION. 

C 

U1=U0*((S1 -LAMDA*S3 ) *COS(BO )+(S2+LAMDA*S4 ) * 

*SIN(B0)+C1 ) 

C 

MUF=XMANG1(1 ./MACHF) 

THTAF=THTAL+XL2-XL 1 
C 

R0=-(  (U  1/SQRT(X2)+SIN(THTAF)  *SIN(MUF)  **2.) 

*/(SIN(2.*MUF) )) 

C 

RADIUS=-1 ./RO 
C 

RETURN 

END 
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SUBROUTINE  CNTLIN(K, S,I) 

C 

C  CALCULATES  CONDITIONS  AT  THE  CENTERLINE  BOUNDARY  CONDITION 
C 

REAL  M(3,45) 

INTEGER  S 

DIMENSION  X(3,45),T(3,45),R(3,45) 

COMMOM/MARRAY/X , M, R,T 
COMMON/LASTV/XML, TL, AL 
COMMON/SHAPE/X7 , CO , R7 , TO , RO , GO , XMO 
C 

WRITE (6, 7) 

7  FORMAT (/T20'*»*  CNTLIN  ) 

C 

X 1 =X ( I , K-S+ 1 ) 

R 1 =R (I ,  K-S+1 ) 

XM1 =M(I ,K-S+1 ) 

T1=T(I , K-S+1 ) 

A1sXMANG2( XMLOC ( XM 1 , GO ) ) 

XM9=1.01 

A4=A1 

T4-T1/2. 

XM4-XM1 

R4=R1/2. 

2760  T9=TAN(A4-T4) 

X3-X 1+R1/T9 
G1=FNG(A4,TL) 

Q1sFNQ(A4,XML) 

C 

C  COMPATIBILITY  EQUATION  FOR  THE  CENTERLINE 
C 

XM3=XM1+(T1-G1)/Q1 

IF(ABS(XM3-XM9).LT.  .0001 )G0  TO  2870 
C 

c 

c 

XM9=XM3 

XM4=(XM1+XM3)/2. 

A3=XMANG2( XMLOC (XM3, GO) ) 

A4=(A1+A3)/2. 

GO  TO  2760 
C 

C  AT  CENTERLINE  THE  RADIAL  C0-0RDINATE=0  AND  THETA-0 
C 

2870  X(I+1, K-S+1 )=X3 
R(I+1 , K-S+1 ) =0 . 0 
M(  1+ 1  ,  K-S+1 )=XM3 
T(I+1 , K-S+1 )=0.0 
C 

RETURN 

END 

C 

C 

C 

c 

c 
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SUBROUTINE  CROSPT(I,J) 

REAL  M(3,45) 

DIMENSION  R(3,45) ,T(3 , 45) ,X(3 , 45) ,XH(5 ) ,XI(5 ) ,XK(5 ) ,XJ(5 ) 
C 

COMMON/MARRAY/X , M, R , T 
COMMON/CNETPT/X 3 , R3 , XM3 , T3 
COMMON/GCROSS/XH , XI , X J , XK 
COMMON/CBKPTS/X 1 , X2 , R 1 , R2 , XM 1 , XM2 , T 1 , T2 
WRITE (6, 7) 

7  FORMAT (/T20'*»*  CROSPT  ***' ) 

C 

C  CHARACTERISTICS  HAVE  CROSSED  INTERPOLATE  TO  DETERMINE 
C  A  NEW  POINT. 

C 

S8=(XH(2)-XH( 1 ) )/(XK(2 )-XK( 1 ) ) 
S9=(XH(4)-XH(3))/(XK(4)-XK(3)) 

S6  =  (XH(3)*-XH(1)+XK(1)«S8-XK(3)*S9)/(S8-S9) 

S7=XH(1 )+S8*(S6-XK(1 ) ) 

C 

WRITE(6,10)S6,S7 

10  FORMAT(/T20,'X(C)=' , F10. 4 ,2X, ' R( C) = ' , FI  0.4) 

C 

V1=XI(1)+(XI(2)-XI(1) )*(36-XK( 1 ) )/(XK(2 )-XK( 1 ) ) 
V2=XI(3)+(XI(4)-XI(3) )®(S6-XK(3) )/(XK(4 )-XK(3 ) ) 

V  3=X J ( 1 )  +  (XJ(2)-XJ( 1 ) )*(S6-XK(1 ) )/(XK(2 )-XK( 1 ) ) 
V4=XJ(3)+(XJ(4)-XJ(3) )®(S6-XK(3) )/(XK(4)-XK(3) ) 

C 

X1=S6 

R1=S7 

XM1=(V 1+V2)/2 . 

T1=(V3+V4)/2. 

IF(XH(2).EQ.O.O)GO  TO  4100 
X2=XK(5 ) 

R2=XH(5 ) 

XM2=XI(5) 

T2=XJ(5) 

C 

CALL  NETPT 
GO  TO  4150 
C 

4100  WRITE(6 , 20) 

20  FORMAT (/T20' NEW  POINT  IS  ON  AXIS') 

C 

C  ITERATE  IF  POINT  ON  CENTERLINE 
C 

R3=0.0 

T3=0.0 

A 1 r XM ANG2 ( XMLOC ( XM 1 , GO ) ) 

XM9  =  1.01 

A4=A1 

T4=T1/2. 

XM4=XM1 

R4=R1/2. 

C 

5540  T9=TAN( A4-T4) 

X3=X1+R1/T9 
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G1=FNG(A4ST4) 

Q1sFNQ(A4sXM4) 

XM3=XM1+(T1-G1)/Q1 
IF(ABS(XM3-XM9).LT..0001)G0  TO  4150 

XM9=XM3 

XM4=(XM1+XM3)/2. 

A3-XMANG2 ( XMLOC ( XM3 , GO ) ) 
A4=(A1+A3)/2. 

GO  TO  5540 


150  X(IS J)sX3 
R(I,J)-R3 
M(I,J)=XM3 
T(I s J)sT3 
T3D-T3S57" 2958 
C 

WRITE (6 , 30)X  3 , R3 , XM3 , T3D 

30  F0RMAT(/T10,X3=,F10.4,2X#,R3='F10.4J2X,'M3=’F10„4> 

»2X,,T3='F10.4) 

RETURN 

END 
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SUBROUTINE  CROSS ( I, K , S, A8) 

INTEGER  S,A8,P,P0,Z0,Z1 
REAL  M(3,45) 

DIMENSION  R(3,45),T(3,45),X(3,45),XI(5), 

*XJ(5) ,XK(5),XH(5) 

COMMON/MARRAY/X , M , R , T 
COMMON/CNETPT/X 3 , R3 , XM3 , T3 
COMMON/CCROSS/XH , XI , X J , XK 
COMMON/CBKPTS/X 1 , X2 , R1 , R2 , XM1 ,XM2,T1 ,T2 
C 

C  CHECKS  IF  CHARACTERISTICS  OF  THE  SAME  FAMILY  HAVE  CROSSED 
C  IF  CHARACTERISTICS  HAVE  CROSSED  INTERPOLATE  TO  DETERMINE  A 
C  A  NEW  POINT. 

C 

C 

WRITE (6, 7) 

7  FORMAT ( /T20 ' ***  CROSS  ***') 

WRITE (6 , 10)1 ,K, S 

10  FORMAT(/T20,,I=,I3,2X'K=,I3,2X'S='I3) 

11=1+1 

IF(I.EQ.  DRETURN 
KS1=K-S+1 

IF(A8.EQ.1)KS1=K-S-1 

C 

DO  20  J=2 , KS1 

C 

J1=J-1 
J2=J  1 
J3=J 
J4-J-1 

IF( A8. NE . 1 )G0  TO  5 
J2=J 
J3=J+1 
J4=J+1 
C 

5  IF( A8 . EQ. 2 . AND. X(I ,J).LT.X(I1,J1))G0  TO  20 
IF(A8. EQ. 1 . AND.X(1 ,J) .LT.X(I1 ,J) )G0  TO  20 
C 

WRITE (6, 30) 

30  FORMAT (/T20* CHARACTERISTICS  HAVE  CROSSED') 

JO=J 

XK( 1 )=X(I-1 , J2) 

XH( 1 ) =R(I- 1 , J2) 

XI(1 )=M(I-1 ,J2) 

XJ(1 )=T(I-1 ,J2) 

C 

XK(2)=X(I,J) 

XH(2 ) =R( I , J ) 

XI(2 ) =M( I , J) 

XJ(2)=T(I,J) 

C 

XK(3 ) =X( I , J1 ) 

XH(3 ) =R(I , J1 ) 

XI(3)=M(I , J 1 ) 

XJ(3)=T(I , J1 ) 

XK(4)=X(I1,J2) 


43 


XH(4 ) =R(I 1 , J2) 

XI(4)=M(I1, J2) 

XJ(4)=T(I1,J2) 

c 

XK(5)-X(I-1,J3) 

■XH(5)=R(I-1,J3') 

XI(5)=M(I-1 s  J  3 ) 

XJ(5)=T(I-1,J3) 

C 

4710  WRITE (6 , 40) (XK(IJ) ,IJ=1 ,5) 

40  FORMAT ( 1X,5F12.4) 

CALL  CROSPTd , J) 

C 

X(I,J)=X3 

R(I,J)=R3. 

M(I,J)=XM3 

T(I,J)=T3 

C 

X1sX3 
R1=R3 
XM1=XM3 
T 1  sT3 

IF(XH(2) . LT.O,Q)GQ  TO  4710 
X2=X(I,J+1) 

R2sR(I,J+1) 

T3-T(I , J+1 ) 

IF(A8,EQ.2)XM2-M(IJJ+1) 

IFU8.EQ. 1)XM3=M(I,J+1) 

CALL  NETPT 
C 

X ( I 1 , J4)sX3 
R(I1,J4)=R3 
M( I 1 ,J4)=XM3 
T(I1, J4)=T3 
C 

WRITE (6 , 40)X3, R3?  XM3, T3 

K=K-1 

KS1=K-S-1 

DO  60  P0=J0+1,KS1 

X(I1,P0-1)=X(I1,P0) 

R ( 1 1 ,  PO-1 )=R(I 1 ,P0) 

M( 1 1 , PO- 1 ) =M( 1 1 , PO) 

T(I 1 ,P0-1 )=T(I 1 ,P0) 

60  CONTINUE 

DO  70  ZlsJO-1 ,K-S 

WRITE (6 , 80)11 , Z 1 ,X(I1,Z1) ,R(I1,Z1) ,M(I1,Z1) ,T(I 1,Z1) 
80  FORMAT (/T30* CORRECTED  POINTS' ,/ ,213, 4F12. 4) ■ 

70  CONTINUE 

20  CONTINUE 

RETURN 
END 
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$PROG  FFREAD 

SUBROUTINE  FFREAD  (HOL, FLO, LH, LF) 

q  * a»» »«**#»#»«*»***»*»*»»»#«*»**#***#*»* a **»»»* 

INTEGER  HOL 

LOGICAL  TRANSP , DIGIT , COMENT , DELIO , DELT05 
DIMENSION  FORM( 6 ) , IBUF(  20) , H0L(20) , FL0(20) 

DATA  NWORDS ,  NIT,  NOT,  NCHARS  /  20,  5,  6,4  / 

TRANSP ( K ) =K . EQ . IBLANK . OR . K . EQ . ICOMMA . OR . K . EQ . I EQUAL 

DIGIT (K)  =K . GE. IZERO. AND. K . LE. ININE 

COMENT ( A ) =IGC( A , 1 ) . EQ. ILETTC . AND. IGC( A , 2 ) . EQ. IBLANK 


PURPOSE 
C 

C  READS  ONE  CARD  IN  FREE  FIELD  FORMAT 

C  THE  $  SIGN  MAY  BE  USED  TO  CANCEL  ITEMS  OR  AN  ENTIRE  LINE 

C  ANY  COMMENT  CARDS  ENCOUNTERED  IN  THE  DATA  DECK  ARE 

C  PRINTED,  SKIPPED  AND  NOT  INTERPRETED 
C 

C  ARGUMENTS  I=INPUT 
C  0=0UTPUT 

C 

C  HOL  —  0  —  WILL  CONTAIN  HOLLERITHS  ENCOUNTERED  ON  CARD 

C  FLO  —  0  —  WILL  CONTAIN  FLOATS  ENCOUNTERED  ON  CARD 

C  INT  —  0  —  WILL  CONTAIN  INTEGERS  ENCOUNTERED  ON  CARD 

C  LH  —  0  —  NUMBER  OF  HOLLERITHS  ON  CARD 

C  LF  --  0  —  NUMBER  OF  FLOATS  ON  CARD 

C  LI  —  0  —  NUMBER  OF  INTEGERS  ON  CARD 

C 

C  DEFINITIONS  —  LIMITATIONS 
C 

C  A  *TRACAR*  (TRANSPARENT  CHARACTER)  IS  A  BLANK,  COMMA  OR  =  SIGN 
C  THE  FOLLOWING- ARE  DELIMITERS 
C  -  BEGIN  OF  CARD 
C  -  END  OF  CARD 

C  -  A  STRING  CONSISTING  OF  ONE  OR  MORE  TRACARS 

C  AN  ITEM  CONSISTS  OF  ONE  OR  MORE  NON-TRACARS  PRECEDED  AND  FOLLOWED  BY  A 
C  DELIMITER. 

C  EACH  ITEM  ON  A  CARD  WILL  BE  INTERPRETED  AS  ONE  OF  THE  FOLLOWING 
C  -  INTEGER 

C  -  FLOAT 

C  -  HOLLERITH 

.  C  ANY  ITEM  STARTING  WITH  +  -  .  OR  DIGIT  WILL  BE  INTERPRETED  AS 
C  -  INTEGER  IF  IT  CONTAINS  NO  . 

C  -  FLOAT  IF  IT  CONTAINS  ONE  . 

C  -  HOLLERITH  IF  IT  CONTAINS  MORE  THAN  ONE  . 

CALL  OTHER  ITEMS  WILL  BE  INTERPRETED  AS  HOLLERITHS 

C  GOOD  INTEGERS  —  1  +5  -208  0034  -051 

C  BAD  INTEGERS  —  ++5  -6+3  2A 

C  GOOD  FLOATS  —  1.E+5  2.-3  -6.5  +006.5  .2 

C  BAD  FLOATS  —  1 . A+5  2.-3X  +6.6—6 

C  HOLLERITHS  —  ..  A  12. .6  +3.-4  XAQ  X5  8  A$A  $$$$X 

C  CANCELED  ITEMS—  A$  43$  12.  .6$  XXXX$  1$  12$  123$  1234$  $$  $$$ 
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C  CANCELED  LINE  —  2.8  4.2  1264  SHIFT  $ 
C 

I PLUS  -IGC(1H+,1) 

ILETTC=IGC( 1 HC ,  1 ) 

IMIN  =IGC(1H-,1) 

I  DOLL  =IGC'(1H$,1) 

I DOT  =IGC( 1 H . , 1 ) 

IBLANK=IGC( 1 H  .1) 

ICOMMA=IGC(1H, ,1) 

IEQUAL=IGC( 1H=, 1) 

IZERO  =IdC( 1HO, 1 ) 

ININE  sIGC( 1H9, 1 ) 

C  INITIALIZE  COUNTS  AND  READ  CARD 

1020  LH=0 
LF=0 
LI-0 

READ  ( 1 , 1060) (IBUF(I) ,1=1 ,NWORDS) 
1060  FORMAT  (20A4) 

IF  (C0MENT(IBUF(1)))  GO  TO  1020 
12=0 


FIND  II  =  FIRST  COLUMN  OF  FIELD 

1100  11=12 
1120  11=11+1 

IF  (I1.GT.80)  RETURN 
IF=IGC(IBUF ,11 ) 

IF  (TRANSP(IF) )  GO  TO  1120 

FIND  12  =  LAST  COLUMN  OF  FIELD 


12=11 

1140  12=12+1 

IF  (I2.GT.80)  GO  TO  1160 
IL=IGC(IBUF ,12) 

IF  ( .NOT. TRANSP(IL) )  GO  TO  1140 
12=12-1 
GO  TO  1180 
1160  12=80 

C  IW  =  FIELD  WIDTH 

1180  NX=I1-1 

IW=I2-I 1+1 
NCH=NX+IW 
DO  1190  1=1,6 
1190  FORM(I)  =4H 

IL=IGC(IBUF ,12) 

IF  ( IL . EQ . IDOLL . AND . I 1 . EQ . 12 )  GO  TO  1020 
IF  (IL.EQ. IDOLL)  GO  TO  1100 
IF=IGC( IBUF  ,11 ) 

IF  ( IF . EQ. I PLUS)  GO  TO  1200 
IF  (IF. EQ. IMIN)  GO  TO  1200 
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IF  (IF.EQ.IDOT)  GO  TO  1200 
IF  (DIGIT(IF) )  GO  TO  1200 
GO  TO  1240 

C  COUNT  THE  DOTS 

1200  ND0TS=0 

DO  1220  1=11,12 

IF  (IGC(IBUF,I) .EQ.IDOT)  ND0TS=ND0TS+1 
1220  CONTINUE 

IF  (NDOTS.EQ.O)  GO  TO  1 440 
IF  (NDOTS.EQ. 1 )  GO  TO  1340 

C  HOLLERITH 

1240  IW=MINO(IW,NCHARS) 

NCH=NX+IW 

IF  (NX.EQ.O)  GO' TO  1280 
ENCODE  (FORM, 1260)NX,IW 
1260  FORMAT  (  1H( ,  12,  4HX, 1A ,  12,  1H)  ) 

GO  TO  1320 

1280  ENCODE  ( FORM, 1 300 )IW 

1300  FORMAT  (  3H(1A,  12,  1H)  ) 

1320  LH=LH+1 

DECODE  ( IBUF ,FORM)HOL(LH) 

GO  TO  1100 


FLOAT 


1340  IF  (NX.EQ.O)  GO  TO  1380 
ENCODE  (FORM, 1 360 )NX,IW 

1360  FORMAT  (  1H(,  12,  3HX,F,  12,  3H.1)  ) 

GO  TO  1420 

1380  ENCODE  (FORM, 1400) IW 

1400  FORMAT  (  2H(F,  12,  3H.1)  ) 

1420  LF=LF+1 

DECODE  (IBUF, FORM) FLO(LF) 

GO  TO  1100 

INTEGER 


1440  IF  (NX.EQ.O)  GO  TO  1 480 
ENCODE  (FORM, 1460 )NX,IW 
1460  FORMAT  (  1H(,  12,  3HX,I,  12*  1H)  ) 

GO  TO  1520 

1480  ENCODE  (FORM, 1500)IW 

1500  FORMAT  (  2H(I,  12,  1H)  ) 

1520  CONTINUE 

DECODE  (IBUF ,FORM)INT 

LF=LF+1 

FLO(LF) =INT 

GO  TO  1100 

END 

C  «»»»«#**»»«**««««***#»**»*»«}>»**#««»**»«*»*#*«»*#** 

FUNCTION  IGC  (ISTR.N) 

C  a*##***###*#**#***#***#******#*#*****’**###*#**#*#*****### 
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DIMENSION  ISTR(I) 

DOUBLE  PRECISION  ISUB 
DATA  ISUB  jNCHARS  /  3HIGC,4  / 

C— 

C  THIS  FUNCTION  RETRIEVES  THE  N-TH  CHARACTER  IN  THE  STRING 

C  <ISTR>  AND  STORES  ITS  OCTAL  6-BIT  CODE  IN  THE  6  LEAST 
C  SIGNIFICANT  BITS  OF  (IGC>. 

C  NOTE  THIS  ROUTINE  CONTAINS  INTERDATA  8/32  DEPENDENT  STATEMENTS 

MWORD  =  Y'OOOOOOFF' 

ICHAR=N- 1 

IWORD-ICHAR/N CHARS 
IRELAsX  CHAR-I WORDsN  CHARS 
IRELAs  (NCHARS-1  )-IRE’LA 
ISHIF=8§IRELA 
JWORD=ISTR(IWORD+1) 

JWORD=ISHFT( JWORD, -ISHIF) 

IGC  =  I AND (JWORD, MWORD) 

RETURN 

END 
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SUBROUTINE  FLOWIN(G AMP , MLP , THTALP , MFP ) 

REAL  MLP, MFP 

DIMENSION  H0L(20) , FL0(20) 

COMMON/RADS/P  I, RAD 
COMMON/VAR/VARd  7) 

THIS  SUBROUTINE  ALLOWS  INPUT  OF  PROTOTYPE 
NOZZLE  FLOW  CONDITIONS 

WRITE( 1 ,10) 

10  FORMAT(/T20' PROTOTYPE  FLOW  SPECIFIED' 

®//T20'ENTER  GAMMA  OF  PROTOTYPE') 

CALL  FFREAD(HOL, FL0,LH,LF) 

GAMP=FL0( 1 ) 

C 

WRITE (1 ,20) 

20  FORMAT (/T20' ENTER  PROTOTYPE  NOZZLE  EXIT  MACH')  ' 

CALL  FFREAD(HOL, FL0,LH,LF) 

MLP=FL0( 1 ) 

C 

WRITE(1 ,30) 

30  FORMAT(/T20' ENTER  PROTOTYPE  NOZZLE  EXIT  ANGLE') 

CALL  FFREAD(HOL,FLO,LH,LF) 

THTALP-FLOO  ) 

C 

WRITEO  ,40) 

40  FORMAT(/T20' ENTER  PROTOTYPE  JET  SURFACE  MACH  NUMBER') 
CALL  FFREAD(HOL, FLO, LH, LF) 

MFP=FLO( 1 ) 

C 

VAR( 1 ) =GAMP 
VAR(2 ) =THTALP 
VAR(3 ) sMLP 
VAR(4 ) =MFP 
C 

THTALB=THTALP/RAD 

C 

RETURN 

END 
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SUBROUTINE  INTERM 

DIMENSION  H0L(20) , FL0(20) ,XI( 11) ,XJ(  1  1 ) ,XK( 11) ,XH( 11) 
COMMON/C 1005/XH , XI , XJ  ,  XK 
1330  WRITE ( 1 s 10) 

10  FORMAT(/T20*ENTER  SEGMENT.*) 

CALL  FFREAD(HOL, FLO,LH,LF) 

IPsFLO( 1 )+1 
WRITE(1 ,20) 

20  F0RMAT(/T20' ENTER  RADIUS') 

CALL  FFREAD(HOL, FLO, LH,LF) 

QsFLO(l) 

U=XJ(IP-1) 

UI-COS(U) 

U2=U1+(Q-XI(IP-1))/XK(IP) 

.  U.2=ATAN(SQRT(1  .-U2SU2)/U2) 

HPKP-XH(IP“1  )+XK(IP){i(SIN(U)-SIN(U2) ) 

U2D=U2»57.3 

WRITE(6 , 30)Q , U2D , HPKP 

30  F0RMAT(.1X,*  R=  *  ?F10.4, '  THETAs  *  ,F10.M  , '  X=',F10.4) 
GOTO  1330 
END 


50 


SUBROUTINE  ITEHAT(GAMM, THTAF , THTALM , RADIUS , MFP , MFM , MLP , MLM) 
R  EAL  L AMD A , MFM , ML , ML2 , ML  3 , ML  P , MLM , MFP 
COMMON/RADS/PI, RAD 
COMMON/VAR/VAR(17) 

C 

C  THIS  SUBROUTINE  CALCULATES  A  VALUE  FOR  THEY A-L-MODEL 
C  BASED  ON  MODEL  GAMMA  AND  A  GUESS  AT  MACH-L-MODEL 
C 

C  ■ 

ML=MLP 

C 

LAMDA=3QRT( (GAMM-1 . ) /(GAMM+1 . ) ) 

Q=THTAF 

RAD2=RADIUS 

X 1 sXMSQR ( ML , GAMM) 

X2=XMSQR (MFM , GAMM) 

CALL  TANGLE ( X 1, XL 1.LAMDA) 

CALL  TANGLE ( X 2 , XL2 , L AMDA ) 

THTALM-Q-XL2+XL 1 
THDEG=THTALMSRAD 
WRITE (6 ,20)ML, THDEG 

20  FORMAT ( /T20 ' ML-M- 1  =  * F 1 0 . -4  ,  *  THETA-L-Ms  ' F 1 0 . 4  ) 

C 

CALL  CALCPL(GAMM, ML , THTALM ,MFM , THTAF , RADIUS ) 

C 

C  SECOND  PASS 
C 

ML2-ML 
RAD3=RADIUS 
ML=MLP+2 . * ( MFM-MFP ) 

C 

X 1 =XMSQR ( ML , GAMM) 

CALL  TANGLE (XI, XL 1.LAMDA) 

THTALM=Q-XL2+XL 1 
THDEG=THTALM#RAD 
WRITE (6 , 40)ML , THDEG 

40  FORMAT( /T20 ' ML-M-2= ' F 1 0 . 4 , ' THETA-L-M-2= ' F 1 0 . 4 ) 

ML3=ML 

C 

CALL  CALCPL(GAMM, ML , THTALM ,MFM , THTAF , RADIUS ) 

C 

C 

C  THIRD  PASS 
C 

50  ML3=ML 

RAD4=RADIUS 

ML=ML+(  RAD4-RAD2) *(ML3-ML2) /(  RAD3-RAD4 ) 

IF(ML. LT. 1 . )ML= (ML2+1 . ) /2 . 

X 1 =XMSQR ( ML , GAMM) 

CALL  TANGLE ( XI, XL 1.LAMDA) 

THTALM=Q-XL2+XL 1 
THDEG=THTALM*RAD 
WRITE(6,60)ML, THDEG 

60  FORMAT ( /T20 ' ML-M= ' F 1 0 . 4 , 2X , ' THETA-L-M= ' F 1 0 . 4 ) 

C 

C  TEST  FOR  CLOSURE 
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IF( ABS(RAD2~RAD4) . LT. .01 )G0  TO  70 

ML2=ML3 

RAD3-RAD4 

CALL  CALCPL(GAMM,ML, THTALM,MFM,THTAF, RADIUS) 
GO  TO  50 

END  CONDITION  MET 


70  MLM=ML 

TRAD~THTAFsRAO 

C 

VAR(5)”TRAD 

VAR(6)=RADIUS 

RETURN 

END 
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SUBROUTINE  JINIT 


C 

C  INITIALIZES  THE  ARRAYS  FOR  THE  METHOD  OF  CHARACTERISTICS 
C  TO  SOLVE  FOR  ACTUAL  PLUME  SHAPE 
C 

REAL  M(21,11),LO,L9,M1,M6,MA(3,45) 

DIMENSION  X(2 1 , 1 1 ), R(2 1 , 1 1 ) ,T(2 1 , 1 1 ) 

* ,XH( 1 1 ) ,XI( 1 1 ) ,XJ( 1 1 ) ,XK(  1 1 ) 
*,RA(3,45),XA(3,45),TA(3,45) 

COMMON/C1 005/XH , XI , X J , XK 
COMMON/INJET/MEXIT,M6,KO,NO 
COMMON/ARRAY/X , R,M,T 
COMMON/SHAPE/X7 , CO , R7 , TO , RO , GO , XMO 
COMMON/INDICE/I , N , K , IS , IP 
COMMON/MARRAY/XA ,MA , RA , TA 
C 

C  SET  INITIAL  CHARACTERISTIC  VALUES  FROM  THE  NORMALIZED 
NOZZLE  EXIT  CONDITIONS. 


DO  100  Jsl ,K0 
X(1 ,J)=XA(1 ,J) 

R( 1 , J)=RA( 1 , J) 

M(  1 , J ) =MA( 1  ,  J) 

T( 1 , J )=TA( 1 , J) 

100  CONTINUE 
C 

LO=SQRT( (GO+1 . ) /(GO-1 . ) ) 

L9=(FNO(FNB(M6,GO) ,LO,GO)-FNO(FNB(XMO,GO) ,L0,G0)+T( 1 , 1 ) ) *57. 295 
WRITE(6,10)L9 

10  FORMAT (///T50' PR ANDTL-MEYER  EXPANSION', 

*//T10'INITIAL  SLOPE  (DEG)  =  ',F10.4) 

1  =  1 
■  IP=2 

XH( 1 )=0. 0 
XI(1)=1. 

XJ(1)=L9*3. 14159/180. 

XK( 1 )=0. 0 
C 

C  SETS  UP  INITIAL  VALUES  IN  ARRAYS  BASED  ON  EXPANSION  AT  LIP 
C  FROM  LIP  EXIT  MACH  NUMBER  TO  JET  SURFACE  MACH  NUMBER. 

C  THE  I  VALUES  ARE  INCREMENTED  FROM  1  TO  NUMBER  OF  NET  POINTS. 

C  ALL  INITIAL  X  VALUES  =  0  AND  ALL  INITIAL  R  VALUES  =1 

C 

WRITE(6 , 40) 

40  F0RMATO HI //,T25' INITIAL  VALUES  FOR  METHOD  OF  CHARACTERISTICS', 

*/T25'TQ  SOLVE  ACTUAL  PLUME  SHAPE') 

J  =  1 

S  TEP= ( M6 -XMO ) / ( FLOAT ( N - 1 ) ) - . 000 1  - 
DO  20  XM=XMO , M6 , STEP 
M1=FNB(XM ,G0) 

.  0=FN0(M1 ,LO,GO) 

0D=0*57 • 3 

WRITEC6 ,30)XM,M1 ,0,0D 

30  F0RMAT(7X,'M=  'F10.4,'  M*=  ' FI 0 . 4 , ’  OMEGA  ='F10.4,'  DEG=  'F10.4) 
X(I , J)=0. 
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R(I,J)=1. 

M( I , J ) =M 1 

T(I, J)=0-FNO(FNB(XMO,GO) ,LO,GO)+T( 1 , 1 ) 
1  =  1+1 

20  ■  CONTINUE 

RETURN 
END 
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SUBROUTINE  KLGLEV(X,Y, J) 

C 

C  THIS  SUBROUTINE  SOLVES  THE  TRANSONIC  FLOW  REGION  IN  THE 
C  NOZZLE  THROAT  SECTION.  CALCULATIONS  ARE  MATE  BEGINNING  AT 
C  THE  THROAT  (X=0,R=1)AND  CONTINUING  UNTIL  THE  MACH  NUMBER  ’ 

C  EXCEEDS  1.025  AT  WHICH  POINT  THE  PROGRAM  SWITCHES  TO 
C  METHOD  OF  CHARACTERISTICS  SOLUTION  TO  SOLVE  FLOW  FIELD  OUT 
C  TO  THE  NOZZLE  EXIT. 

C 

REAL*8  L1,L2,L3,L4",L5,N1,N2,N3,N4,N5,N6,Z0,Z7,Z8,Z9 
DIMENSION  B(45),C(45) 

COMMON/ CKLEV/B , C 

C0MM0N/SHAPE/X7 , CO , R7 » TO  *  RO , GO , XMO 
WRITE(6 ,7 ) 

7  FORMAT (/T20'*»*  KLGLEV  ***') 

Z0=X*SQRT(2.*C0/(G0+1 .) ) 

L1=Y*Y/2.-. 25+ZO 

L2=((2.*G0+9.)*Y*»4-(4.*G0+15.)*Y*Y)/24. 

L2=L2+(10.*GO+57.)/288.+ZO*(Y*Y-.625)-(2.*GO-3.)*ZO*ZO/6. 

L3= (556. *G0**2+1 899. *GO+3231.)*Y**6/1 0368.- (388. *GO*GO 
&+ 1 233 • *G0+1 953 • ) *Y**4 
L3-L3/2304. 

L3=L3+( 304. *G0*G0+858.*G0+1 269. )*Y*Y/1 728.- (2708. *GO*GO+ 
&7839.*G0+14211 .) 

L3=L3/82944. 

L4=(52.*GO*GO+99.*GO+375.)*Y**4/384.-(52.»GO*GO+99.*GO+303.) 
&*Y*Y/192.+(92.*G0*G0+l80.*G0+639.)/1 152. 

L5= ( ( 1 3. *G0-27 . ) /48. - (5 . *G0-5 . ) *Y*Y/8 . ) *ZO*ZO+ZO*ZO“ZO* 
*((4.*G0*G0-57.*G0+27.)/144.) 

L3=L3+Z0*L4+L5 

B(J)=1 . +L 1/( CO+1 . ) +(L 1+L2) /( CO+1 . ) **2+(L1+2 . *L2+L3) / 
&(C0+1.)**3 
N1=Y*Y*Y/4.-Y/4.+Y*Z0 

N2=(8.»G0+15.)*Y**5/72.-(20.*G0+45.)*Y**3/96.+ 

&(28.*G0+75.)*Y/288. 

N2=N2+Z0*((4.*G0+9.)*Y**3/12.-(4.*G0+9.)*Y/12.) 

Z8=Y**5/1 3824 . 

Z7=6836.*GO*GO 

N3= (Z7+1 6695. *G0+1 42 11.)*Y**7/82944. -(3380. *GO*GO+ 

&8703 .*G0+7875 . ) *Z8 

N3=N3+(3424 .*G0*G0+91 83 . *G0+8964 . ) *Y**3/1 3824 . 

N3=N3-(7 100 .*G0*G0+1 9575 .*G0+20745 . ) *Y/82944 . 
N4=(556.*GO*GO+1113.*GO+981.)*Y**5/1728.- 
&( 388 . *G0*G0+801 . *G0+693 . ) *Y**3/576 . 

N4=N4+(304 .*G0*G0+645 .*G0+549 • ) *Y/864 . 
N5-Z0*Z0*((52.*G0*G0+3.*G0-33.)*I**3/192.~ 
&(52.*G0*G0+27.*G0-9.)*Y/192.) 

N5=Z0**3*(G0+1 . ) *Y/4 . 

N3=N3+Z0*N4+N5-N6 

Z9=SQRT( (GO+1 . ) /(2 .* (CO+1 . ) ) ) 

Z8=1.5*N1+N2 

C( J ) =Z9* (N 1 /( CO+1 .)+(Z8)/( CO+1.) **2+( 15. *N1/8.+ 

&2.5*N2+N3) / ( CO+1 . ) **3 ) 

RETURN 

END 
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SUBROUTINE  USTPL(THTAF ,  RADIUS) 


C 

C 

C  THIS  SUBROUTINE  PRINTS  THE  PLUME  SHAPE 

C 

C 

WRITE(8 ,15) 

WRITE (9 , 10) 

10  FORMAT(  1  HI  ,//,40X,' PLUME  SHAPE  CALCULATED  FROM  -JOHANNESEN ' 
*' -MEYER  THEORY* 

»/T56,'X* ,T66,'R' ,T76, ’THETA') 

15  FORMAT( 1H1 ,// ,5X, ’PLUME  SHAPE  CALCULATED  FROM  JOHANNESEN* 
a* -MEYER  THEORY* 

*/T21,,X',T31,,R',T4l, 'THETA') 

RAD2=Q„0 

DO  30  S=0.  ,5.  ,  .1 
X=SIN(THTAF)~ (S/RADIUS) 

IsXMANGI (X) 

S1-RADIUSg(SIN(THTAF)--SIN(T) ) 

R=1 .+RADIUS*(COS(T)-COS(THTAF) ) 

TDEG=Tf57- 3 
WRITE(8 , 25)  S1sRfTDEG 
WRITE(9,20)  SI , R,TDEG 
20  FORMAT(T51 »3F10.4) 

25  FORMAT (T 1 6, 3F1 0.4) 

IF(RAD2, GT. R) RETURN 
RAD2=R 

30  CONTINUE 

END 
C 


56 


o  o 


SUBROUTINE  LISTPR(MF, GAMMA) 

REAL  LAMDA.MF 
COMMON/VAR/VAR(17) 

THIS  SUBROUTINE  CALCULATES  AND  PRINTS  THE  PRESSURE  RATIO 
C 

LAMDA-SQRT( (GAMMA-1 . ) /(GAMMA+1 . ) ) 

X=XMSQR(MF, GAMMA) 

Ps ( 1 . -XSLAMDASLAMDA) 89 (GAMMA/ ( GAMMA- 1 . ) ) 

VAR(7)=P 

RETURN 

END 
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SUBROUTINE  MODIN (G AMM, MFM , GAMP , MFP , IERR ) 

REAL  MFM 

DIMENSION  HOL(20) , FL0(20) 

COMMON/VAR/VAR(17) 

C 

C  THIS  SUBROUTINE  ALLOWS  INPUT  OF  MODEL  GAMMA 
C  AND  MODEL  FINAL  MACH  NUMBER 
C 

WRITE (1 , 10) 

10  FORMAT(/T20' ENTER  GAMMA  OF  MODEL') 

C ALL  FFREAD ( H OL , FLO , LH , LF ) 

GAMMsFLO( 1 ) 

C 

WRITE (1 ,20) 

20  FQRMAT(/T20» ENTER  ;  1 -INPUT  FINAL  MACH  OF  MODEL' 

*/T2 9* 2 -CALCULATE  WITH  WEAK  SHOCK  RELATIONS' 

*/T29* 3 -CALCULATE  WITH  STRONG  SHOCK  RELATIONS') 
CALL  FFREAD(HOL, FLO, LH, LF) 

I0PT=FL0(1 ) 

GO  10(30,40,50) ,IOPT 
30  WRITE (1,22) 

22  FORMAT(/T20' ENTER  FINAL  MACH  NUMBER  FOR  MODEL') 

VAR( 1 6) =MFM 
RETURN 
C 

40  CALL  WEAKFM( GAMP, MFP ,GAMM, MFM, IERR) 

VAR( 1 6) =MFM 
RETURN 
C 

50  CALL  STRGFM( GAMP, MFP, GAMM, MFM, IERR) 

VAR( 1 6) =MFM 

RETURN 

END 


58 


o  o  o  o  o 


SUBROUTINE  NETPT 
C 

C  THIS  SUBROUTINE  RETURNS  THE  INTERMEDIATE  MOC  NET 
C  VALUES  (X3,R3,XM3,T3)  BASED  ON  THE  1  AND  2  VALUES 
C 

COMMON/ CBKPTS/X 1 , X2,R1 , R2 , XM 1 , XM2 , T1 , T2 
COMMON/CNETPT/X  3 , R3 , XM3 , T3 
COMMON/LASTV/XM5 , T5 , A5 
COMMON/SHAPE/X7 , CO , R7 , TO , RO , GO , XMO 
C 

WRITE (6 ,7) 

7  FORMAT ( /T20' *»*  NETPT  «««' ) 

XM9-1 .01 
C 

THE  A  VALUES  ARE  MACH  ANGLES 

A1=XMANG2(XML0C(XM1 ,G0)) 

A2=XMANG2(XML0C(XM2,G0) ) 

A4-A1 
A5=A2 
T4=T  1 
T5=T2 
XM4-XM1 
XM5-XM2 
R4=R1 
R5=R2 

X3  AND  R3  ARE  THE  POSITION  OF  THE  INTERMEDIATE  NET  POINT 

1530  X3=(R2-R1+X1“TAN(T4-A4)-X2*TAN(T5+A5)) 
&/(TAN(T4-A4)-TAN(T5+A5) ) 

R3-R1+(X3-X1 )*TAN(T4-A4) 

C 

Q1=FNQ(A4,XM4) 

G1=FNG(A4,T4) 

C 

Q2=FNQ(A5,XM5) 

C 

C  CHECK  IF  POINT  IS  ON  CENTER  LINE 
C 

IF( R2 . EQ. 0 . 0)G0  TO  2630 
F2=FNF(A5,T5) 

C 

C  XM3  IS  THE  MACH  NUMBER  OF  THE  INTERMEDIATE  NET  POINT 
C  T3  IS  THE  THETA  OF  THE  INTERMEDIATE  NET  POINT 
C 

XM3=((T1-T2)+G1*(R3-R1)/R4+F2S(R3-R2) 

&/R5+Q 1SXM1+Q2SXM2) /(Q 1+Q2) 
T7=Q2*(XM3-XM2)-F2*(R3-R2)/R5 
1680  T3=T2+T7 

C 

C  CHECK  IF  LOCAL  MACH  NUMBER  CLOSE  TO  SONIC 
C 

IF(ABS(XM9-XM3) .LT. .0001 )RETURN 

XM9=XM3 

T4=(T1+T3)/2. 
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T5=(T2+T3)/2. 

XM4=(XM1+XM3)/2. 

XM5  =  ( XM2+XM3 ) /2 „ 

A3=XMANG2(XML0C(XM3,G0) ) 

A4=(A1+A3)/2. 

A5=(A2+A3)/2, 

R5=(R3+R2)/2, 

RMR3+R1 )/2. 

GO  TO  1530 

POINT  IS  ON  CENTERLINE  CALUCULATE  MACH  AND  THETA 
C 

2630  XM3=  (Q2tXM2/2  .+T1-s-Q1sXM1+G  1a(R3-R1 )  /R4 )/ 

& ( Q  2/2 . +Q 1 ) 

T7=Q2*(XM3~XM2)/2. 

GO  TO  1680 
END 
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SUBROUTINE  NOZEXT(I ,MEXIT) 

REAL  M(3 ,45) 

DIMENSION  X(3,45) ,R(3,45) ,T(3,45) 

COMMON/MARRAY/X , M , R , T 
COMMON/C1 820/X8, R8, XM8 , T8 
COMMON/SHAPE/X7, CO , R7 , TO , RO , GO , XMO 
COMMON/VAR/VA'R(17) 

C 

C  AFTER  NOZZLE  EXIT  CONDITION  PASSED',  MODIFY.  CHARACTERISTICS 
C  TO  LIP  OF  NOZZLE. 

C 

W1=0.0 

10  IF( MEXIT . EQ.O) WO =(RO-R (1-1, 1))/(R(I+1,1)-R (1-1,1)) 

IF(MEXIT.EQ. 1)W0=(XM0-XML0C(M(I-1 , 1) ,G0))/ 

*(XML0C(M(I+1 ,1) ,G0)-XML0C(M(I-1 ,1) ,G0)) 

C 

C  LINEARLY  INTERPOLATE  EXIT  CONDITIONS 
C 

X(I , 1 ) =X (I- 1 , 1)+(X(I, 1)-X(I-1 , 1) )*wo 
R(I,  1  )=R(I-1 ,1)  +  (R(I, 1)-R(I-1 , 1) )*W0 
M( I , 1 ) =M( I- 1 , 1 ) + (M( I , 1 ) -M( I- 1 , 1 ) ) *W0 
T(I,1)=T(I-1,1)+(T(I,1)-T (1-1,1))* WO 
C 

CALL  WALL(I) 

W0SW1=ABS(W0-W1 ) 

C 

WRITE (6 , 769 )V/0 ,  W1 ,  W0SW1 
769  FORMAT ( 5  X , 3F 1 5 . 5 ) 

IF( ABS(WO-WI) .LT.  .01)G0  TO  20 
X(I+1 , 1)=X8 
R(I+1 , 1 )=R8 
M(I+1 , 1)=XM8 
T(I+1 , 1 )=T8 
W1=W0 
GOTO  10 
C 

20  XM88=XMLOC(XM8 ,G0)  . 

T8D=T8*57. 3 
VAR( 1 3)=X8 
VAR( 1 4) =R8 
C 

RETURN 

END 
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SUBROUTINE  PLULIS 
COMMON/VAR/VAR( 1 7 ) 


C 

C 

C 

10 


20 


30 


40 


50 


THIS  SUBROUTINE  OUTPUTS  THE  PROGRAM  INPUTS  AND  RESULTS 


WRITE(9 , 10) (VAR(I) ,1=1 ,7 ) 

FORMAT ( 1  HI , T56 , ' MODEL  NOZZLE  DESIGN ' 

1//T4 5' PROTOTYPE' 

2/  T45 '  — - * 

3/T42' SPECIFIC  HEAT  RATIO, ........................ 

4 /T4 2' NOZZLE  ANGLE.  ............................... 


5/T42'EXIT  MACH  NUMBER. ........................... 

6/T42' JET  SURFACE  MACH  NUMBER. . ................... 

7/142' INITIAL  SLOPE  OF  JET  PLUME. 


8/T42!INITIAL  RADIUS  OF  CURVATURE  OF  JET  PLUME.... 

9/T42'PRESSURE  RATIO. ............................. 

WRITE(9,20)(VAR(I) ,1=8, 16) 


FORMAT  ( //, T45 ,' MODEL ',/, T45 
1/T42' SPECIFIC  HEAT  RATIO. ........ . 

2/T42 ' THROAT  RADIUS  OF  CURVATURE... 


3/T42' BEGINNING  AXIAL  LOCATION  OF  CONICAL  SECTION. 
4/T42' BEGINNING  RADIAL  LOCATION  OF  CONICAL  SECTION 
5/T42!NOZZLE  ANGLE ,  DEG.  .........................  . 


6/T42'NOZZLE  LENGTH. ................. 

7/T42* NOZZLE  EXIT  RADIUS. ............ 

8/T42'EXIT  MACH  NUMBER, .............. 

9/T42' JET  SURFACE  MACH  NUMBER ....... . 

IF( VAR( 1 7) . LE.  2.)GO  TO  50 
WRITE (8 , 30) (VAR(I ) , 1= 1 ,7) 

FORMAT ( 1  HI , T21 , 'MODEL  NOZZLE  DESIGN' 
1//T1 O' PROTOTYPE' 


2/  T10' - ' 


3/T7 'SPECIFIC  HEAT  RATIO .......................... 

4/T7 'NOZZLE  ANGLE. ................................ 

5/17 ' EXIT  MACH  NUMBER ............................. 

6/T7 ' JET  SURFACE  MACH  NUMBER. ..................... 

7/17 'INITIAL  SLOPE  OF  JET  PLUME,..,.,.,........... 

8/T7 'INITIAL  RADIUS  OF  CURVATURE  OF  JET  PLUME.,... 
9/T7 'PRESSURE  RATIO. .............................. 

WRITE (8 ,40)(VAR(I) ,1=8,16) 

F0RMAT(//,T10, 'MODEL'  ,/,T10,'— '  , 

1/T7 'SPECIFIC  HEAT  RATIO. ......................... 

2/T7 'THROAT  RADIUS  OF  CURVATURE..,,...,,......,... 

3/T7 'BEGINNING  AXIAL  LOCATION  OF  CONICAL  SECTION.. 
4/T7 'BEGINNING  RADIAL  LOCATION  OF  CONICAL  SECTION. 


5/T7 'NOZZLE  ANGLE,  DEG....... 

6/T? 'NOZZLE  LENGTH.  .......... 

7/T7 'NOZZLE  EXIT  RADIUS...... 

8/T7 ' EXIT  MACH  NUMBER ........ 

9/T7 'JET  SURFACE  MACH  NUMBER. 
VAR(5)=VAR(5)/57.3 
CALL  LISTPL( VAR(5 ) , VAR(6 ) ) 
CALL  WRITPL 
RETURN 
END 


. ' F8. 3 , 
. s  F8 . 3  j 
.'F8.3, 
• ' F8. 3, 
. ' F8. 3 , 
.'F8.3, 
. »F8.3) 


.'F8.3, 
. ' F8 . 3 , 
.’F8.3, 
. '  F8 . 3 , 
. ' FQ . 3 , 
• ' F8. 3 , 
=  *  F8 . 3 , 
. ' F8. 3 , 

.’F8.3) 


' F8 . 3 , 
1 F8 . 3 , 
' F8. 3, 
' F8. 3, 
*  F8 . 3 , 
'F8.3, 
' F8.3) 


*  F8 . 3 , 
’F8.3, 

*  F  8 . 3 , 
'F8.3, 
'  F8 . 3 , 
'F8. 3, 
'F8.3, 
'F8.3, 
'F8.3) 
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SUBROUTINE  PLUMIN(THTAFP , RADP) 

DIMENSION  H0L(20) ,FLO(20) 

COMMON/RADS/PI, RAD 
COMMON/VAR/VAR0  7) 

C 

C  THIS  SUBROUTINE  ALLOWS  INPUT  OF  PLUME  SHAPE 
C 

C  THTAFP-THETA-F  OF  PROTOTYPE 
C  RADP-INITIAL  RADIUS  OF  PROTOTYPE 
C 

WRITE( 1 , 10) 

10  F0RMAT(/T20 'ENTER  THETA-F  AND  INITIAL  RADIUS' 

*/T22'0F  PROTOTYPE') 

C 

CALL  FFREAD(HOL, FLO,LH,LF) 

THTAFP-FL0O  ) 

RADP-FLO(2 ) 

IF( LF , EQ=  2)  GO  TO  20 
CALL  FFREAD(HOL, FLO, LH, LF) 

RADP=FL0(1) 

C 

20  VAR(5 ) =THTAFP 

VAR(6 ) =RADP 
C 

THTAFPsTHTAFP/RAD 

C 

RETURN 

END 
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SUBROUTINE  PLUMPT 

REAL  M(21, 1 1 ) , Ml , M2 , M3 , M5 , M6 

DIMENSION  X(2 1 , 1 1 ) , R(2 1 , 1 1 ) ,T(2 1 , 1 1 ) ,XH( 1 1 ) ,XK( 1 1 ) 

* ,XI( 1 1 ) ,XJ( 1 1) 

COMMON/I N JET/MEXIT , Mo,KO,NO 
COMMON/ ARRAY/X,R,.M,T 
COMMON/ C 1 005/XH , XI , X J , XK 
COMMON/INDICE/I , N , K , IS , IP 
COMMON/SHAPE/X7 , CO , R7 , TO , RO , GO , XMO 
J=1 

X2sX(N , J+1 ) 

R2=R(N,J+1) 

M2=M(N,J+1) 

T2=T(N , J+1 ) 

A2sXMANG2(XML0C(M2,G0)  ) 

A3=XMANG2(XMLOC(M(N,J) ,GO) ) 

R1=R(Nf J) 

X1=X(N,J) 

M1=M(N,J) 

T1=T(N,J) 

A5sA2 

M3=FNB(M6,GO) 

A3=XMANG2(XMLOC(M3,GO) ) 

T5-T2 

M5=M2 

R5sR2 

T9=T1 

X3=X2 

X3=(R1  -R2-i-X2sTAN(T5+A5 )-X  1aTAN(T9 )  )/(TAN(T5+A5)-TAN(T9 ) ) 
1105  R3=R2+(X3~X2)*TAN(T5+A5) 

Q2=1 ./ (M5*TAN(A5) ) 

F2  =S  IN  ( A5 )  *S  IN  ( T5 )  /SIN  ( A-5  +T5 ) 

T3=T2+Q2»(M1-M2)-F2S(R3-R2)/R5 

T9=(T1+T3)/2. 

X7  -  (R 1  -R2+X28T  AN(T5+A5  )-X  1'3TAN(T9)  )/(TAN(T5+A5)-TAN(T9 ) ) 

IF(ABS(X3-X7)  .LT.  ,001 )G0  TO  1170 

X3=X7 

T5=(T2+T3)/2. 

M5“(M2+M3)/2. 

A5-(A2+A3)/2. 

R5-(R2+R3)/2. 

GOTO  1105 

1170  RO=(R1-R3)/(COS(T1)-COS(T3)) 

C 

T3D=T3957.3 

X(N+1 , 1 ) =X3 

XH(IP)=X3 

XK(IP)=RO 

R(N+1,1)=R3 

XI(IP)=R3 

M(N+1 , 1)=M3 

T(N+1,1)=T3 

XJ(IP)=T3 

RETURN 

END 
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SUBROUTINE  PRTOPT 
DIMENSION  HOL(20) , FLO (20) 

COMMON/VAR/VAR07) 

C 

C  SETS  UP  FILES  AND  LOGICAL  UNITS  FOR  OUTPUT 
C 

1  WRITE( 1,10) 

10  FORMAT (/T20' ENTER  OUTPUT  PREFERENCE'/, 

*/T20'PRINTER. .................  1 '  , 

9/T20'VARIAN. .................. 2 '  , 

8 /T20' CONSOLE. . ................  3', 

*/T20' CONSOLE  AND  PRINTER . 4'  , 

a/T20' CONSOLE  AND  VARIAN. _ 5') 

C 

CALL  FFREAD(HOL, FLO,LH,LF) 

IPRT0P=FL0(1) 

VAR(17)=IPRT0P 

IF(IPRT0P.NE. 1  .AND.  IPRTOP. NE, 4)G0  TO  20 
CALL  CLOSE(9  > ISTAT) 

CALL  OPENW(9, 'PR:', 2,0,0, ISTAT) 

IF(IPRT0P.EQ.1)G0  TO  100 
C 

20  IF( IPRTOP. NE. 2  .AND.  IPRTOP. NE. 5 )G0  TO  30 

CALL  CLOSE (9 i ISTAT) 

CALL  QPENW(9, 'PLOT: 1 ,2,0,0, ISTAT) 

IF( ISTAT.NE. 0)GO  TO  40 
IF( IPRTOP. EQ. 5 )G0  TO  30 
GO  TO  100 
40  WRITE (1,22) 

22  FORMAT (/T 15' 99  ERROR  WITH  VARIAN  -  CHECK  AND  RE-ENTER88 ’ ) 

GO  TO  1 
C 

30  CALL  CLOSE (8, ISTAT) 

CALL  0PENW(8 , 1  CON : 1 ,4, 0,0, ISTAT) 

C 

100  RETURN 
END 
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SUBROUTINE  PRWALL(I) 


THIS  SUBROUTINE  PRINTS  OUT  CONDITIONS  AT  POINTS 
C  ALONG  THE  WALL( PRINTS  TO  UNIT  6) 

C 

REAL  M(3S45) 

DIMENSION  R(3»45) ,T(3,45) »X(3,45) 
COMMON/MARRAY/X , M  s R » T 
COMMON/SHAPE/X7 , CO , R7 , TO , RO , GO „ XMO 
C 

WRITE (6,7) 

7  FORMAT ( /T20’  PRWALL  ii3a') 

D 1=XMANG2(XML0C(M(I-1 , 1 ) ,GO) ) 
D2=XMANG2(XML0C(M(I, 1) ,GO) ) 
D3=“COS(D1)+CQS(D2)»M(I , 1 )/M(I™1 , 1) 

D4sSQRT(  (X(I-1,1)-X(I,1))i‘#2  + 
»(R(I-1,1)-R(I,1))**2)/R(I-1,1) 

C 

U1=D3*M(I-1 , 1)/D4 
D5=D3/D4+SIN(T0)»(SIN(D2))i»»2 
D6=D5/(2,9(COS(D2) )**2) 

C 

RORsR(I- 1,1) 

XOR=X(I-1 , 1) 

c 

WRITE(6 , 10)U 1 ,D6, ROR,XOR 
10  FQRMAT( /T1Q, 'U1=’,F10.4,5X,' ACCEL™ ’ , FI  0. 4 , 
*/T10'AT  W ALL  POINT  WHERE  R/R*  =',F10.4, 

8 /T1 O' AND  X/R  = ' , F10. 4) 

C 

RETURN 

END 
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SUBROUTINE  RESET (K,S) 


C 

C  THIS  SUBROUTINE  STORES  THE  VALUES  OF  THE  CHARACTERISTIC 
C  POINTS  JUST  CALCULATED  AND  CURRENTLY  IN  THE  1=3  POSITION 
C  INTO  THE  1=1  POSITION.  THIS  SUBROUTINE  CAL.,', ED  EVERY  OTHER  . 
C  PASS.  THEREFORE,  TWO  PLANES  OF  CHARACTERISTICS  ARE  CARRIED 
C  ALONG  ALWAYS  AND  THE  THIRD  PLANE  ONLY  TEMPORARILY. 

C 

REAL  M(3,45) 

INTEGER  S 

DIMENSION  X(3,45) ,R(3,45) ,T(3,45) 

COMMON/MARRAY/X , M, R ,  T 
WRITE (6, 7) 

7  F0RMAT(/T20'®»*  RESET  *»»• ) 

C 

DO  10  J=1 ,K-S+1 
X(1 ,J)=X(3,J) 

R(1 ,J)=R(3,J) 

M(1,J)=M(3,J) 

T(1 ,J)=T(3,J) 

10  CONTINUE 

C 

RETURN 

END 
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SUB  ROUTINE  SETEXT (P , A8 , 1 , Q , K , S , N , IFINI ) 

C 

c 

REAL  M(3sM5) 

INTEGER  S,P,A8,Q 

DIMENSION  D(3,45) ,E(3,45) ,F(3,45) ,G(3,45) 

»sX(3s45) sR(3,45) 

COMMON/MARRAY/X ,M,R,T 
COMMON/C2480/D ,  E, F, G 
C 

C  SETS  UP  CHARACTERISTICS  AFTER  EXIT  REACHED 
C 

C  NOTE:  P<0  AFTER  EXIT,  QsO  UNTIL  EXIT  THEN  INCREASED 
C  BY  1  EVERY  OTHER  PASS 
WRITE (6, 7) 

7  FQRMAT(  /T20 1  •**  SETEXT  **■*')■ 

C 

IF(A8.EQ.2.  .AND.P.LT.O.)  Q=Q+1 
11=1+1 

WRITE(6, 10)11 

10  FORMAT( /T20 ! ROW  (1+1)=', 14) 

DO  20  J=1 ,K-S 
IF(J„LT.Q)GQ  TO  20 
IF(J.NE.Q,OR.P,GE,0)G0  TO  2540 
N=N+ 1 

D ( 1 ,N)=X(I+1,J) 

E(1,N)=R(I+1,J) 

F(  1 ,N)=M(I+1,J) 

G(  1  ,N)=T(I+1  ,J) 

C 

C 

2540  TDEG=T(I 1 , J) *57. 3 

WRITE(6,30)J,X(I1 , J) »  R( 1 1 , J) ,M(I1 , J) ,TDEG 
30  FORMAT(1X,I4,4F10.3,2X, 'SETEXT1 ) 

IF((2*(K-S)-1).EQ.N)G0  TO  40 
20  CONTINUE 

RETURN 

40  CALL  STOREX(K,S) 

IFINI=1 

RETURN 

END 

C 

C 

C 

C 
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SUBROUTINE  SETPTS(K , S , I) 

THIS  SUBROUTINE  SETS  UP  THE  TWO  POINTS  NEEDED  TO 
CALCULATE  THE  NEXT  INTERMEDIATE  NET  POINT 

INTEGER  S 
REAL  M(3,45) 

DIMENSION  X(3,45),T(3,45),R(3,M5) 

COMMON/MARRAY/X , M, R , T 

COMMON/ CBKPTS/X 1 , X2 , R 1 , R2 , XM 1 , XM2 , T 1 , T2 

COMMON/CNETPT/X 3 , R3 , XM3 , T3 

WRITE (6, 7) 

FORMAT(/T20' SETPTS  *#*• ) 


DO  10  J= 1 , K-S 
X1-X(I,J) 
X2=X(I,J+1) 
R1=R(I,J) 
R2=R(I, J+1) 
XM1 =M( I , J) 
XM2=M(I, J+1) 

T  1  =T( I , J ) 
I2=T( I , J+1 ) 
CALL  NETPT 
X(I+1 ,J)=X3 
R(I+1, J)=R3 
M(I+1 , J ) =XM3 
T(I+1, J)=T3 
0  CONTINUE 


RETURN 

END 
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SUBROUTINE  SETPT2 
REAL  M(21,11),M1,M2,M3 
DIMENSION  X(2  1 , 1 1 ) ,  R(2  1 , 1 1 )  ,  T.(2  1 , 1 1 ) 
COMMON/ CBKPTS/X 1 , X2 , R1 , R2  ,  Ml , M2 , T1 , T2 
COMMON/ CNETPT/X  3 , R3 , M3  ,  T3 
CQMMON/ARRAY/X , R , M,  T 
COMMON/INDICE/ I , N , K , IS, IP 
DO  10  J=1,K-1 
X1=X(I+1,J) 

X2=X(I , J+1 ) 

R1=R(I+1,J) 

R2=R(I ,J+1 ) 

Ml =M( 1+ 1 , J ) 

M2=M( I , J+1 ) 

TUT  (1+1  ,J) 

T2=T(I,J+1) 

CALL  NETPT 
X(I+1 , J+1 )=X3 
R(I+1,J+1)=R3 
M(I+1 , J+1)=M3 
T(I+1 , J+1 )=T3 
10  CONTINUE 
RETURN 
END 
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SUBROUTINE  SETPT3 

REAL  M(2 1,11), Ml, M2, M3 

DIMENSION  X(2 1 , 1 1 ) , R(2 1 , 1 1 ) ,T(2 1 , 1 1 ) 

COMMON/CNETPT/X 3 , R3 , M3 , T3 

COMMON/CBKPTS/X 1 , X2 , R1 , R2 , Ml , M2 , T1 , T2 

COMMON/ARRAY/X , R,M,T 

COMMON/INDICE/I , N , K , IS , IP 

DO  10  J=1,K-2 

X1=X(N+1 ,J) 

X2=X(N,J+2) 

R 1 =R(N+ 1 , J) 

R2=R(N , J+2) 

M1=M(N+1 , J) 

M2sM(  N  ,  J$-2 ) 

T1=T(N+1 ,J) 

T2=T(N , J+2) 

IF( M2 . LT .  1 , ) RETURN 
CALL  NETPT 
X(N+i,J+1)=X3 
R(N+1 ,  J-s-1  )=R3 
M(N+1 , J+1 )~M3 
T(N+1 , J+1 )=T3 
10  CONTINUE 
RETURN 
END 
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SUBROUTINE  STOREX(K,S) 

C  .  ■ 

C  CALLED  AFTER  SUFFICIENT  CHARACTERISTICS  PAST  LIP 
C  ARE  CALCULATED.  THE  EXIT  CHARACTERISTICS  ARE 
C  NORMALIZED  AND  STORED. 

C 

INTEGER  S 
REAL  M(3,45) 

DIMENSION  X(3 , 45) ,1(3,45) ,R(3, 45) ,0(3,45) ,E(3,45) ,F(3,45) 

9 ,G(3 , 45) 

C 

C0MM0N/C2480/D , E, F,G 
COMMON/MARRAY/X ,M,R,T 
COMMON/SHAPE/X 7 , CO , R7 , TO , RO , GO , XMO 
C 

WRITE(6,7) 

7  FORMAT(/T20!  STOREX  «**') 

WRITE (6, 10) 

10  FORMAT (/T20* NOZZLE  EXIT  CHARACTERISTIC  NORMALIZED  AND  STORED*) 
WRITE (6 ,20) 

20'  F0RMAT(/T23'N'T30'X(1 ,N) *T40'R(1 ,N) * T50 *M( 1 , N) * T60 1 T( 1 , N)  RAD* ) 

C 

KS2=2*(K-S)-1 
DO  30  N  =  1 ,KS2 

X(1 ,N)=(D(1 ,N)-D(1 ,1))/E(1  ,1) 

R(1 ,N)=E(1 ,N)/E(1 ,1) 

T(1,N)=G(1,N) 

M(  1  ,  N)  sF(  1  ,  N) 

C 

WRITE(6,40)N,X(1 ,N) , R(1 , N) , M( 1 ,N) , T( 1 , N) 

C 

30  CONTINUE 

C 

£XITM=XMLOC( M( 1,1), GO ) 

WRITE(6 , 60) EXITM 

60  FORMAT ( / 1 0X , ' EXIT  MACH  NUMBER  AT  LIP  =  *,F12.4) 

40  FORMAT(I24,4F10.3) 

C 

RETURN 

END 
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SUBROUTINE  STRGFM( G AMP , MFP , GAMM , MFM , IERR ) 
REAL  MFP, MFM 


C 

C  SOLVES  FOR  THE  FINAL  MACH  NUMBER  OF  THE  MODEL ,  ALSO 
C  REFERRED  TO  AS  THE  JET  SURFACE  MACH  NUMBER. 

C  THE  STRONG  SHOCK  APPROXIMATION  FOR  MATCHING  THE  SUPERSONIC 
INVISCID  STREAMLINE  DEFLECTION-PRESSURE  RISE  RELATION  IS  USED. 

A=(2.«GAMP»MFP*MFP-GAMP+1 . ) /(GAMP+1.) 

B=Aa (GAMM+1 . ) +GAMM-1 . 

C=B/ (2 . *GAMM) 

IF(C.GT.O.O)GO  TO  10 
WRITE(1 ,20) GAMP, GAMM, MFP 
20  F0RMAT(//T10' A  SOLUTION  FOR  THE  JET  SURFACE  MACH  NUMBER' 

*/T10'DOES  NOT  EXIST  FOR  THE  SPECIFIED  CONDITIONS*8 ' 
*/T15’GAMMA  PROTOTYPE  -  'FI 0.2, 

«/T15'GAMMA  MODEL  '=  'F10.2, 

S/T15’MACH  MODEL  =  'F10.2) 

IERR-1 

RETURN 

C 

10  MFM=SQRT(C) 

C 

RETURN 

END 
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SUBROUTINE  WALL(I) 

C 

C  THIS  SUBROUTINE  SOLVES  FOR  THE  WALL  BOUNDARY 
C 

REAL  M( 3 , 45) 

DIMENSION  X(3,45),R(3,45),T(3,45) 
COMMON/MARRAY/X , M, R , T 
COMMON/SHAPE/X7 , CO , R7 , TO , RO , GO , XMO 
COMMON/L ASTV/XML , TL , AL 
COMMON/C1 820/X8 , R8 , XM8 , T8 
C 

WRITE (6 ,7) 

7  FORMAT (/T20' WALL  ***•) 

C 

X2=X(I , 1 ) 

R2=R(I,1) 

XM2=M(I,1) 

T2=T(I , 1 ) 

A2  =XM ANG2 ( XMLOC ( XM2 , GO ) ) 

XM9=1 .01 

A5sA2 

T5=T2 

XM5=XM2 

R5-R2 

1920  T9=TAN(T5+A5) 

C 

C  CHECK  IS  CONICAL  SECTION  REACHED 
C 

IF(X2„LE,X7.AND.X3.LE,X7)GO  TO  2060 
C 

C  CALCULATE  X,R, THETA  IN  CONICAL  SECTION 
C 

2590  X3=(R2-R7+X7STAN(T0)-X2ST9)/(TAN(T0)-T9) 
R3=R7+(X3-X7)STAN(T0) 

T3-T0 

GO  TO  2121 
C 

C  CALCULATE  X,R, THETA  IN  CIRCULAR  THROAT  SECTION 
C 

2060  X0=C0+1.-R2+T9*X2 

X9=T9aX0/(1 .+T9ST9) 

X3=X9-SQRT(X9sX9-s-(C0*C0-X0*X0)/(T9*T9+1.)) 
IF(X3.GT.X7)GO  TO  2590 
R3=1 .+C0»(1 .-SQRT(1 .-X3aX3/(CO»CO))) 
T3=ATAN( 1 ./SQRT( C0*C0/(X3*X3)~1 . ) ) 
T3D=T3S57. 2958 
WRITE (6, 10)X3,R3»T3D 
10  FORMAT(3F12.5 ,2X, ' WALL' ) 

2121  Q2=FNQ(A5,XML) 

F2=FNF(A5,TL) 

XM3-XM2+(T3»T2+F28(R3-R2)/R5)/Q2 
IF(ABS(XM9-XM3) .LT.  .0001 )G0  TO  2220 
XM9-XM3 
T5=(T2+T3)/2. 

XM5  =  ( XM2+XM3 ) /2 . 

A3=XMANG2(XMLOC(XM3,GO)) 
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A5=(A2+A3)/2. 
R5-(R3+R2)/2. 
GO  TO  1920 
C 

2220  X8=X3 

R8-R3 
XM8-XM3 
T8=T3 
RETURN 
END 


SUB  ROUTINE  WEAKFM( GAMP , MFP , GAMM , MFM , IERR) 

REAL  MFP 8 MFM 
C 

C  SOLVES  FOR  THE  FINAL  MACH  NUMBER  OF  THE  MODEL,  ALSO 
C  REFERRED  TO  AS  THE  JET  SURFACE  MACH  NUMBER , 

C  THE  WEAK  SHOCK  APPROXIMATION  FOR  MATCHING  THE  SUPERSONIC 
C  INVISCID  STREAMLINE  DEFLECTION-PRESSURE  RISE  RELATION  IS  USED. 
C 

A=GAMP/GAMM*MFPsMFP/SQRT(MFPsMFP-1 . ) 

IF( A ,GS .2 . )GQ  TO  10 
WRITE(1 ,20)GAMP, GAMM, MFP 

20  FORMAT ( //T 10 ’ 88  A  SOLUTION  FOR  THE  JET  SURFACE  MACH  NUMBER’ 

*/T10'D0ES  NOT  EXIST  FOR  THE  SPECIFIED  CONDITIONS8*' 
a/T 15' GAMMA  PROTOTYPE  =  'F10.2, 

S/T15’GAMMA  MODEL  =  ’F10.2, 

*/T15’MACH  MODEL  =  ’F10.2) 

IERRsI 

RETURN 

C 

10  MFM=SQRT((AsA+A*SQRT(A*A-4.) )/2.) 

C 

RETURN 

END 
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SUBROUTINE  WRITPL 

DIMENSION  XK( 1 1 ) , XI ( 1 1) ,XJ(1 1) ,XH(1 1) 

COMMON/INDICE/I , N , K , IS , IP 
COMMON/C  1 005/XH , XI , X J , XK 
C 

WRITE  OUT  RESULTS 

WRITE(8, 15) 

WRITE (9, 10)  • 

10  FORMAT(1H1 ,//T38'MODEL  PLUME  SHAPE  CALCULATED  BY  METHOD  OF' 
*'  CHARACTERISTICS' 

*/T56,'X' , T66 , ' R ' ,T76, 'THETA' ) 

15  FORMAT ( 1  HI ,//T2 'MODEL  PLUME  SHAPE  CALCULATED  BY  METHOD  OF' 
«'  CHARACTERISTICS' 

9/T21,'X' ,T31 ,'R' ,T4l , 'THETA' ) • 

DO  30  1=1, IP- 1 

WRITE(8,25)XH(I) ,XI(I ) , XJ(I ) «57 . 3 
WRITE  (9 ,20)XK(I),XH(I),XI(I)-,XJ(I)st57.3 
20  F0RMAT(T51 .3F10.4) 

25  FORMAT(T16,3F10.4) 

30  CONTINUE 

C 

RETURN 

END 
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SUBROUTINE  WRITRW 


C 

C  CALLED  TO  WRITE  OUT  INTERMEDIATE  VALUES  WRITES  GO  TO 
C  LU  6  ONLY  CALLED  FOR  DEBUGGING 
C 

REAL  M( 21, 11) 

DIMENSION  X(21, 1 1 ) , R(2 1 , 1 1 ) ,T(2 1 , 1 1 ) 

COMMON/ ARRAY/X ,  R ,  M,  T 
COMMON/INDICE/I ,N,K, IS, IP 
C 

I 1 “ 1+ 1 

WRITE (6, 10)11 

10  FORMAT( /T20 1  ROW  '=  *,I3) 

C 

DO  20  J s 1 , K 

TD=T(I1,J)»57-.3 

WRITE (6  »  30) J ,X(1 1 s J ) , R(1 1 , J) , M( 1 1 , J) , TD 

30  FORMAT (IX, 13, 4F 12. 4) 

20  CONTINUE 

C 

RETURN 

END 
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SUBROUTINE  YINCR(P , MEXIT, I , K , S) 

REAL  M(3,45) 

INTEGER  S,P 

DIMENSION  X(3,45) ,R(3,45) ,T(3,45) 

COMMON/MARRAY/X , M, R,T 
COMMON/C  1 820/X8 , Rff, XM8 , T8 
COMMON/SHAPE/X7 , CO , R7 , TO , RO , GO , XMO 
C 

C  MOVES  THE  CURRENTLY  CALCULATED  CHARACTERISTIC  POINTS  UP  ONE 
C  IN  THE  Y  POSITION  OF  THE  ARRAYS  (J) 

C 

WRITE (6, 7) 

7  FORMAK /T20'  YINCR  ****  ) 

KS1=K-S+1 

DO  10  J=KS1,1,-.1 

X(I+1,J+1)=X(I+1,J) 

■R(I+1,J+1)=R(I+1,J) 

M(I+1,J+1)=M(I+1,J) 

T(I+1 , J+1 )=T(I+1 , J) 

10  CONTINUE 

C 

C  SET  WALL  CONDITIONS  INTO  J=1  POSITION 
C 

X(I+1 , 1)=X8 
R(I+1,1)=R8 
M(I+1,1)=XM8 
T(I+1,1)=T8 
C 

C  CHECK  FOR  CONICAL  SECTION 
C 

K=K+2 

IF(X(I+1,1).LT.X7)G0  TO  2410 
CALL  PRWALL(I) 

C 

C  P-1  BEFORE  NOZZLE  EXIT  AND  P=-1  AFTER  NOZZLE  EXIT 
C 

2410  IF( P . LT . 0 ) R  ETURN 
C 

C  TEST  IF  SPECIFIED  EXIT  MACH  NUMBER  EXCEEDED 
C 

IF ( MEXIT. NE. 1 )GOTO  40 
IF(XMLOC(M(I+1,1),GO).GT.XMO)GO  TO  30 
RETURN 
C 

C  TEST  IF  SPECIFIED  EXIT  RADIUS  EXCEEDED 
C 

40  IF(R(I+ 1 , 1 ) . LE. RO) RETURN 

C 

30  Ps-1 

CALL  NOZEXT ( I , MEXIT ) 

RETURN 

END 
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FUNCTION  XMANG1 (XMACH) 


C 

C  THIS  FUNCTION  RETURNS  THE  MACH  ANGLE 
(SAME  AS  SIN  ALPHAsl/M) 

THE  INPUT  XMACH  IS  1/M 

XMANG1 =ATAN (XMACH/SQRT( 1 .-XMACh*XMACH+1 , £-50) ) 
RETURN 
END 


FUNCTION  XMSQR( XMACH, GAMMA) 

THIS  FUNCTION  CALCULATES  THE  VALUE  OF  M»  SQUARED 
WHERE  M®-V/Cs-V/V». 

NOTE:  THIS  ADIABATIC  RELATIONSHIP  CAN  BE  FOUND 
IN  CHAP.  4  PAGE  81  OF  SHAPIRO'S  COMPRESSIBLE  FLUID  FLOW 

XMSQFU  (  ( (GAMMA-s-1  . )  /2  . )  ®XMACH^XMACH)  / ( 1  .  +  ( (GAMMA-1 . )  sXMACHaXMACH) 
&/2.) 

RETURN 
END 
C 

C 

SUBROUTINE  TANGLE(XMACH , THETA , LAMDA ) 

REAL  LAMDA 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  PRANDTL-MAYER  FUNCTION 
C 

C  CHAP.  15  SHAPIRO 
C 

XK=  ( ( 1  .-XMACH) /(XMACH-1 . /LAMDA/ LAMDA)  )’s3(  .5) 
THETA=ATAN(XK)/LAMDA-ATAN(XK/LAMDA) 

RETURN 

END 
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c 

Q&  §  #  &  $  §  H  $  \k  If  #  ¥i  &  &  #  #  ®  &  M  $  §  &  fc  §  $  $  $  -1  &  H  &  it  $  ft  IHHI  §  #  &  It  #  II  ft  >.HH?  $  It  #  #  &  it  &  ft  ft  it  If  it  #  if 

c 

FUNCTION  XMANG2(X) 

THIS  FUNCTION  RETURNS  THE  MACH  ANGLE 
THE  INPUT  IS  MACH 

XMANG2=ATAN( 1 . /SQRT(X»X-1 . ) ) 

.  RETURN 
•  END 

•;;•  a  *  «  a  s 9  a  »  »  s  •:*<  a » •»  «  »  a ;; » :■!  •»  s  «  a  a «  *  *  »  a  » ,•<  :>  a  s  •»  a  « a  2  « ::•  •;:•  •»  •;;• » ;<  »  »  -;s  a  •»  :;■  a :» s 


C  THIS  FUNCTION  RETURNS  LOCAL  MACH  NUMBER 
C  THE  INPUT  IS  M* 

C 

FUNCTION  XMLOC(X.GO) 

XMLOC=SQRT(2 . ) *X/SQRT(G0+1 (GO-1 . ) *XaX ) 

RETURN 

END 

C 

C»a«as-»  s  »»a«a  aaaaaaa  aijaisasMaasasiaassssaaa&ajiasiaasaasaaMaaaaa 

C 

FUNCTION  FNQ( A ,XM) 

FNQ-1 . /XM/TAN( A) 

RETURN 

END 
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c 

c 

FUNCTION  FNF( A ,T) 

FNFr  SIN ( T ) *S IN ( A ) /SIN ( T+A ) 

RETURN 

END 

C 

C#  »®  mss  8B»»iBi»s5  mb  #*«#«»  »*»M«ss#ii3»ss&»a!H4»«»s##»ii»»s*sa»iH»#»is 
C 

FUNCTION  FNG( A ,T) 

FNG=SIN(T)*SIN(A)/SIN(T-A) 

RETURN 

END 

C 

FUNCTION  FNB(X,GO) 

FNB=XsSQRT( (GO+1 . ) /(2 .+(G0-1 . ) «X*X) ) 

RETURN 
END 
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c 

FUNCTION  FNC(X ,G0) 

FNC-SQRT(  (XSX-1 . )  /( (GO+1 . )  /(GO-1  .  )-X*X.) ) 

RETURN 

END 

C 

c 

FUNCTION  FNO(X , XLO ,G0) 

FNO=XLO»ATAN(FNC(X ,G0) )-ATAN(XLO«FNC(X,GO) ) 

RETURN 

END 
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